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1.0  INTRODUCTION 


Supported  by  AFRLA^A  under  a  three-year  STTR  contractual  effort,  ZONA  Technology  has 
successfully  completed  the  seamless-integration  of  the  ZAERO  module  with  ASTROS 
(Automated  Structural  Optimization  Systems).  ZONA  Technology  (ZONA)  and  its  team 
members:  University  of  Oklahoma,  Universal  Analytics,  Inc.  and  Technion,  have  further 
validated  and  tested  the  ZAERO/ASTROS  (called  ASTROS*).  Utilizing  ZAERO,  ZONA  and 
Technion  have  developed  an  Aeroservoelastic  (ASE)  module  for  ASTROS*  during  the 
contractual  period  (called  ASTROServo). 

ZAERO  is  zona’s  unsteady/steady  aerodynamic  program  that  contains  four  essential  modules, 
namely:  the  AGM  (Aerodynamic  Geometry  Module),  the  3D  Spline  Module,  the  UAIC  (Unified 
Aerodynamic  Influence  Coefficient  Module)  and  the  Flutter  Module.  Central  in  ZAERO  is  the 
UAIC  module,  which  renders  ZAERO  applicable  to  complex  aircraft  configurations  covering  all 
Mach  numbers  ranging  from  subsonic,  transonic,  supersonic  to  hypersonic  flight  regimes. 

The  functionality  of  ZAERO  is  to  provide  ASTROS  a  much  improved  aerodynamic  module 
which,  as  opposed  to  the  existing  DLM/CPM  and  USSAERO  codes,  unifies  all  Mach  number 
and  generates  high-fidelity  wing-body  configurated  aerodynamics  for  advanced  aeroelastic 
analysis/design  applications  in  an  ASTROS/MDO  environment.  Thus,  when  interfaced  with  the 
ASE  module,  ZAERO/ASTROS  or  ASTROS*  can  perform  aircraft  design/analysis  with 
additional  aeroservoelastic  constraints. 

In  this  section,  we  will  briefly  describe: 

•  Overview 

•  ZAERO/UAIC:  A  Unified  AIC-Based  Aerodynamic  Module 

•  ZAERO/UAIC  and  ASE  Modules  in  ASTROS 

•  Other  new  ZAERO  modules  in  ASTROS:  AGM,  3D  Spline  and  Flutter 

1.1  Overview 

For  modem  aeroelastic  methodology,  preference  to  the  Computational  Fluid  Dynamics  (CFD) 
methods  or  the  Aerodynamic  Influence  Coefficient  (AIC)  methods  for  aeroelastic  applications 
has  been  the  subject  of  much  discussion.  Here,  we  refer  to  both  methodologies  as  a  part  of 
computational  aeroelasticity.  In  general,  our  concept  of  computational  aeroelasticity  consists  of 
Aeroelastic  Modeling  Methodology,  which  includes  AIC  methods,  stractural  FEM,  spline 
methods  and  flutter  solution  methods,  etc.,  and  Aeroelastic  Simulation  Methodology,  which 
includes  CFD  methods,  closely-coupled  CFD/CSD  interfacing  method,  etc.  (Figs  1.1,  1.2).  In 
our  estimation,  there  should  exist  little  conflict  in  the  choice  of  these  two  methodologies  if  we 
are  directed  towards  a  holistic  approach  for  aeroelastic  applications.  Rather,  they  should 
compliment  each  other  if  their  practices  could  follow  the  proposed  global  strategy  as  shown  in 
Figs  1.1  and  1.2. 
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Aeroelastic  Modelin; 


Fig  1.2  Computational  Aeroelasticity  for  MDO  Applications 


While  both  methodologies  can  cope  with  configurations  of  complex  geometry  in  high  fidelity, 
AIC  methods  can  provide  expedient  amplitude-perturbed  eigen  solutions  in  the  k-domain  and  the 
s-domain.  Rapid  flutter  solutions  can  be  readily  obtained  by  several  standard  flutter  methods. 
For  this  reason,  AIC  methods  are  preferred  by  industries. 

Moreover,  application  of  AIC  methods  to  aeroservoelasticity  and  the  Multidisciplinary  Design 
Optimization  (MDO)  environment  is  straightforward.  In  terms  of  aeroelastic  applications,  they 
should  provide  selected  critical  conditions  for  CFD  methods  to  fine  tune  the  unsteady 
aerodynamics  in  a  confined  flow  regime,  thus  saving  substantial  computing  effort  in  search  of 
potential  flutter  solutions. 

On  the  other  hand,  the  utilization  of  CFD  methods  is  to  link  up  with  a  structural  FEM  via  a 
closely-coupled  CFD/CSD  interfacing,  such  as  the  BEM  solver  (Ref  1.5),  as  indicated  in  Figs  1.1 
and  1.2.  Clearly,  CFD  methods  are  required  when  more  accurate  solutions  become  mandatory  in 
a  flow  regime  where  nonlinearity  dominates  (e.g.  thick  wing  in  supercritical  flow,  high-angle-of- 
attack  flow  with  vortex  dynamics). 

For  classical-flutter  predictions,  the  flow  nonlinearity  could  be  linearized  through  a  robust 
indicial  method  routine  in  conjxmction  with  a  proposed  modal  AIC  method.  In  this  way,  CFD 
solutions  could  be  carried  over  to  the  k-domain  for  its  subsequent  participation  to 
aeroservoelasticity  and  MDO  applications.  For  applications  in  static  aeroelasticity,  the  proposed 
modal  AIC  method  can  be  an  expedient  means  in  utilizing  CFD  solutions  to  generate  a  flexibility 
correction  to  the  measured  rigid  load.  Therefore,  it  is  clear  that  an  effective  AIC  method  with 
stridently  high-fidelity  modeling  capability  remains  to  be  the  backbone  of  computational 
aeroelasticity. 

Towards  this  end,  a  Unified  Aerodynamic  Influence  Coefficient  (UAIC)  method  for  arbitrary 
wing-body  configurations  has  been  developed  that  covers  the  complete  fli^t  regime  of  subsonic, 
transonic,  supersonic  and  hypersonic  Mach  numbers.  This  imsteady/steady  UAIC  methodology 
has  been  further  established  as  a  computer  module,  known  as  the  ZAERO/UAIC  module.  The 
ZAERO/UAIC  module  is  a  stand-alone  aerodynainic  module,  which  can  be  interfaced  with 
existing  FEM  programs  such  as  NASTRAN  and  ASTROS.  Under  a  two-year  AFRLAVright  Lab 
contractual  support,  a  seamless  integration  of  the  ZAERO/UAIC  module  into  ASTROS  is 
successfully  completed.  Fig  1.3  shows  the  integrated  ASTROS/ZAERO  program  architecture. 
This  theoretical  manual  attempts  to  describe  the  theoretical  formulations  of  ZAERO,  with  special 
emphasis  on  UAIC  formulation  and  its  applications  in  each  flow  regime. 

1.2  ZAERO/UAIC;  A  Unified  AIC  Based  Aerodynamic  Module 

The  ZAERO/UAIC  module  consists  of  four  major  unsteady  aerodynamics  codes  that  jointly 
cover  the  complete  domain  of  all  Mach  number  ranges,  namely  ZONA7U  (formerly  ZONA51U), 
ZONA6,  ZONA7  and  ZTAJC.  As  can  be  seen  in  Fig  1.4,  the  aero  modules  currently  integrated 
within  MSC/NASTRAN  and  ASTROS  only  have  the  purely  subsonic  and  supersonic 
capabilities. 
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The  development  of  the  ZAERO  module  has  been  the  major  endeavor  of  ZONA  Technology  in 
the  last  decade.  The  following  is  a  brief  account  for  the  capability  of  the  computer  codes  within 
the  ZAERO  module. 


1.3  ASTROS/ZAERO  Program  Architecture 


Mach  Range  - >■ 


Fig  1.4  ZAEROAJAIC  Module 
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ZONA6/ZONA7:  Generates  Unsteady  Subsonic/Supersonic  Aerodynamics  for  Aircraft 

Confisurations  with  External  Stores 

Prior  to  1990,  all  unsteady  aerodynamic  methods  for  aeroelastic  computations  were  based  on 
low-order  lifting-surface  models  (e.g.  DLM,  Refs  1.6, 1.7).  The  low-order  formulation  of  these 
lifting-surface  methods  renders  input  complexity  in  dealing  with  general  planform 
configurations.  Most  importantly,  the  aerodynamic  effects  due  to  the  presence  of  bodies  and  due 
to  wing-body  interference  were  largely  ignored.  Meanwhile,  the  coupled  external-store  wing 
flutter,  a  problem  that  is  of  frequent  concern  to  modem  aeroelasticians,  can  no  longer  be  resolved 
by  lifting  surface  modeling  alone.  For  this  reason,  development  of  methods  such  as  ZONA6  and 
ZONA7  are  mandatory  (Refs  1.8, 1.9). 

The  main  features  of  ZONA6  and  ZONA7  include; 

•  ZONA6  generates  steady/unsteady  subsonic  aerodynamics  for  wing-body/aircraft 
configurations  vrith  external  stores/nacelles  including  the  body-wake  effect. 

•  ZONA6  is  based  on  a  higher-order  panel  formulation  for  lifting  surfaces  than  the  Doublet 
Lattice  Method  (DLM).  Cases  studied  confirm  the  ZONA6  robustness  over  the  DLM. 

•  Panel  formulation  of  ZONA7  for  lifting  surface  is  identical  to  that  of  ZONA51  -  now  the 
industrial  standard  method  for  supersonic  flutter  analysis  in  MSC/NASTRAN. 

•  ZONA6  and  ZONA7  can  model  any  combinations  of  planar/nonplanar  lifting  surfaces  with 
arbitrary  bodies  such  as  fuselage,  stores,  tip  missiles,  or  their  combinations. 

•  High-order  paneling  of  ZONA6  and  ZONA7  allows  high-fidelity  modeling  of  complex 
aircraft  with  arbitrary  stores/tip  missile  arrangement. 

For  further  details  in  the  background,  application  examples  and  theoretical  formulation  of 
ZONA6  and  ZONA7,  the  reader  is  referred  to  Section  3.0. 

ZTAIC:  Generates  Unsteady  Transonic  Aerodynamics  for  Liftins  Surface  Systems 

Since  1985,  ZONA  has  been  following  up  on  the  development  of  the  Transonic  Strip  (TES) 
Method  (Refs  1 . 1 0, 1 . 1 1 )  for  unsteady  flow  computations  of  arbitrary  wing  planforms.  The  TES 
method  consists  of  two  consecutive  steps,  to  a  given  nonlinear  Transonic  Small  Disturbance 
Code  such  as  ZTRAN,  namely  the  chordwise  mean  flow  correction  and  the  spanwise  phase 
correction.  The  spanwise  correction  procedure  is  further  enhanced  by  a  recently  uncovered 
generalized  relation  obtained  according  to  Oyibos'  separability  principle. 

Based  on  the  TES  concept,  ZONA's  Transonic  Aerodynamic  Influence  Coefficient  (ZTAIC) 
method  is  developed  to  fully  automate  the  computation  procedure  resulting  in  a  modal-based 
AIC  matrix  (Ref  1.12).  The  computation  procedure  requires  direct  pressure  input  from  a  set  of 
computed  or  measured  data.  Otherwise,  it  does  not  require  airfoil  shape  or  grid  generation  for  a 
given  planform.  Meanwhile,  all  the  mean-flow  shock  jumps  are  properly  included  in  the 
resulting  unsteady  aerodynamics  through  the  AIC  formulation.  The  unsteady  pressures  can 
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readily  be  solved  on  the  surfaces  of  a  lifting  surface  system  according  to  the  following  modal- 
based  AIC  formulation,  i.e. 

{ACp}  =  [MAIC]  {h) 

where  [MAIC]  =  [  ACp ]  ,  <!>  is  the  base-line  modes.  ACp^  is  the  computed 

pressure  due  to  and  h  is  the  given  modes  which  expressed  in  terms  of 

The  main  features  of  ZTAIC  include: 

•  ZTAIC  generates  unsteady  transonic  modal  AIC  using  externally-provided  steady  mean 
pressure. 

•  ZTAIC  adopts  steady  pressure  input  (provided  by  measurement  or  CFD),  whereby: 

-  no  grid  generation  is  required,  and 

-  the  correct  unsteady  shock  strength  and  position  are  ensured. 

•  The  modal  AIC  of  ZATIC  as  an  archival  data  entity  allows: 

-  repetitive  aeroelastic  analysis  and  structure  design. 

-  the  ease  of  application  of  the  K  /  P-K  /  g  methods  for  flutter  analysis. 

•  ZTAIC  is  readily  integrated  with  ZONA6  as  a  unified  subsonic/transonic  AIC  method  for 
complex  aircraft  configurations. 

•  Additional  input  to  ZONA6  amounts  to  only  the  provided  steady  pressure  data. 

For  further  details  in  the  background,  application  examples  and  theoretical  formulation  of 
ZTAIC,  the  reader  is  referred  to  Section  4.0. 

ZONA7U:  Generates  Unified  Unsteady  Hypersonic/Supersonic  Aero-dvnamics  for  Liftins. 
Surface  Systems  and  Wim-bodv  Confizurations 

A  Unified  Supersonic/Hypersonic  Lifting  Surface  Method  has  been  developed  recently  (Refs 
1.13,  1.14).  This  method  combines  the  Supersonic  Lifting  Surface  Theory  (such  as  ZONA51, 
Ref  1.15)  with  a  nonlinear  thickness  correction  matrix  Ey,  based  on  a  composite  third-order 
theory,  which  is  rendered  uniformly  valid  throughout  the  hypersonic/supersonic  regime,  i.e. 

{ACp}  =  [D  +  /iEr  M 

where  D  is  the  linear  supersonic  downwash  matrix  provided  by  ZONA51  and  /r  is  a  switching 
function  that  operates  on  the  nonlinear  thickness  matrix  E  for  compression  and  expansion  waves. 
This  correction  matrix  takes  the  flow  nonlinearity  as  well  as  the  flow  rotationality  due  to  shock 
waves  into  account,  which  covers  both  the  Mach-wave  and  Newtonian  limits.  For  aeroelastic 
applications,  ZONA51U  has  been  applied  to  various  wing  planforms  with  thickness 
distributions.  Superseding  ZONA51U,  ZONA7U  integrates  ZONA51U  into  ZONA7  in  that  the 
lifting  surfaces  are  subject  to  unified  hypersonic/supersonic  aerodynamics. 

The  main  features  ofZONA7U  include: 
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•  Z0NA7U  generates  unified  hypersonic  and  supersonic  steady/unsteady  aerodynamics  for 
wing-body/aircraft  configurations  with  external  stores/nacelles. 

•  Nonlinear  thickness  effects  of  ZONA7U  yields  good  agreement  with  Euler  solution  and  test 
data. 

•  Steady  solutions  approach  linear  and  Newtonial  limits. 

•  Confirms  hypersonic  Mach  independent  principle. 

•  Results/formulation  are  superior  to  Unsteady  Linear  Theory  and  Piston  Theory. 

•  ZONA7U  usually  results  in  more  conservative  flutter  boundaries  than  other  methods. 

•  Unified  with  ZONA7  and  is  therefore  applicable  to  all  Mach  numbers  >1.0. 

•  Additional  input  to  ZONA7  amounts  to  only  wing  root  to  wing  tip  sectional  thickness 
profiles. 

For  further  details  in  the  background,  application  examples  and  theoretical  formulation  of 
ZONA7U,  the  reader  is  referred  to  Section  5.0. 

1.3  ZAERO/UAIC  and  ASE  Modules  in  ASTROS 

According  to  the  ASTROS/ZAERO  program  architecture  (Fig  1.3),  database  entities  (such  as 
MAIC)  generated  by  the  ZAERO  module  are  computed  during  the  ASTROS  preface  phase  and 
need  not  be  recomputed  in  the  ASTROS  analysis/optimization  loop.  Meanwhile,  computation  of 
the  ZAERO  module  is  triggered  by  the  new  bulk  data  entry  MKAEROZ  which  specifies  the 
Mach  number,  reduced  frequencies,  method  flags  and  the  mean  flow  conditions. 

The  development  of  an  aeroservoelasticity  (ASE)  module,  based  on  the  ZAERO/UAIC 
aerodynamics,  and  its  integration  with  ASTROS  has  been  completed.  The  ASE  module  will 
facilitate  the  inclusion  of  multi-input  multi-output  (MIMO)  control  system  effects  on  the 
dynamic  stability  and  response  in  the  ASTROS  multidisciplinary  analysis  and  design 
optimization  software  package.  Its  overall  capabilities  include; 

•  Provide  closed-loop  robust  stability  analysis. 

•  Add  continuous  gust  response  capabilities. 

•  Allow  the  inclusion  of  stability  and  gust-response  constraints  in  structural  design 

optimization.  .  ,  ,  •  v 

•  Allow  the  inclusion  of  user-defined  control  parameters  of  a  given  control  law  m  the 
multidisciplinary  optimization  process. 

•  Export  an  efficient  state-space  representation  of  the  aeroservoelastic  system  for  subsequent 
analysis  and  control  synthesis  with  commercially  available  tools  such  as  MATLAB  and 
MATRIX  X. 

ASE  mandates  the  s-domain  aerodynamics  as  a  base,  which  can  be  obtained  ffoni  the  k-domain 
aerodynamics  by  means  ojf  several  rational  approximation  methods.  Among  existing  metiiods 
the  minimum-state  approach  (MIST)  is  selected  here  because  it  offers  significant  savings  in  the 
number  of  added  states  with  little  or  no  penalty  in  the  accuracy  of  modeling  the  aerodynamic 
forces.  The  minimum-state  approach  converts  the  generalized  aerodynamic  forces  Q{ia))  to  0(5) 
in  the  following  form; 


Q(s')  =  Ao  +  Ai 


sb 

y 


+  A2 


+  D  F'\s)  E 


UJ 
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where  F(s)  =  —I  -  R  and  Ai  are  the  real-value  approximation  matrices,  i?  is  a  diagonal 

matrix  with  distinct  negative  terms  representing  the  aerodynamic  lags,  and  D  and  E  are 
aerodynamic  coupling  matrices  between  the  modal  and  aerodynamic  states. 

For  further  details  in  the  background,  application  examples  and  theoretical  formulation  of  the 
ASE  module  utilizing  ZAERO/UAIC,  the  reader  is  referred  to  Ref  1 . 1 6, 1 . 1 7  and  1.18. 

1.4  Other  New  ZAERO  Modules  in  ASTROS:  AGM,  3D  Spline  and  Flutter 

ZAERO  also  includes  three  other  new  modules  in  ASTROS  (Fig  1.4).  These  are:  the 
Aerodynamic  Geometry  Module  (AGM),  the  3D  Spline  Module  and  the  Flutter  Module.  The 
essential  features  of  these  modules  are  briefly  described  as  follows. 

Aerodynamic  Geometry  Module  (AGM) 

The  AGM  module  is  capable  of  modeling  any  full  aircraft  configuration  with  stores  and/or 
nacelles.  A  complex  aircraft  configuration  can  be  represented  by  the  AGM  module  by  means  of 
wing-like  and  body-like  definitions.  Any  modifications  to  the  AGM  module,  such  as  input 
geometry  enhancements,  will  have  minimaJ  impact  on  other  general  modules. 


Wing-LIke  Components  Body^Like  Components 


Wing-Like 

Components 

Include; 

-  wings,  tails,  pylons, 
launchers  store  fins, 
etc. 


Body-Like 

Components 

Include: 

-  fuselage,  underwing 
stores,  missile  bodies, 
etc. 


Fig  1.5  The  Aerodynamic  Geometry  Module  (AGM)  of  ZAERO 
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3D  Spline  Module 

The  3D  Spline  module  establishes  the  displacement/force  transferal  between  the  structural  Finite 
Element  Method  (FEM)  model  and  the  ZAERO  aerodynamic  model.  It  consists  of  four  spline 
methods  that  jointly  assemble  a  spline  matrix.  These  four  spline  methods  include:  (a)  Thin  Plate 
Spline;  (b)  Infinite  Plate  Spline;  (c)  Beam  Spline  and  (d)  Rigid  Body  Attachment  meftods.  The 
spline  matrix  provides  the  x,  y  and  z  displacements  and  slopes  in  three  dimensions  at  all 
aerodynamic  grids. 


Fig  1.6  The  3D  Spline  Module  of  ZAERO 
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Flutter  Module 


The  ZAERO  flutter  module  contains  two  flutter  solution  techniques:  the  AT-method  and  the  P-K 
method.  The  AT-method  is  a  new  contribution  of  ZAERO  to  ASTROS.  For  further  details  on  the 
background  and  the  theoretical  formulation  of  the  3D  Spline  Module  and  the  Flutter  Module,  the 
reader  is  referred  to  Sections  6.0  and  7.0,  respectively. 

Lastly,  we  remark  that  in  parallel  to  seamlessly  integrated  ZAERO  and  ASE  modules  in 
ASTROS,  the  stand-alone  ZAERO  and  ASE  modules  are  being  developed  which  should  be 
readily  interfaced  with  other  systems  such  as  NASTRAN. 
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2.0  AEROELASTICITY  FOUNDATION  OF  ZAERO 


In  this  section,  we  will  discuss: 

•  fundamentals  of  aeroelasticity. 

•  aeroelastic  matrix  equations  of  ZAERO  for  flutter  analysis. 

•  three  disciplines  that  are  required  to  generate  the  aeroelastic  matrix  equations,  namely  the 
structural  finite  element  method,  the  unsteady  aerodynamic  methods,  and  spline  methods. 

Since  the  structural  finite  element  method  (FEM)  is  a  well-established  methodology,  no  detailed 
description  of  the  formulation  for  the  generation  of  the  structural  matrices  will  be  given.  For  the 
unsteady  aerodynamic  methods,  the  concepts  of  amplitude  linearization,  modal  approach  and 
frequency-domain  formulation  for  the  generation  of  Aerod)mamic  Influence  Coefficient  (AIC), 
leading  to  an  eigenvalue  problem  for  the  flutter  analysis,  will  be  discussed  in  detail. 

2.1  Fundamentals  of  Aeroelasticity 

Aeroelastic  response  of  flight  vehicle  is  a  result  of  the  mutual  interaction  of  inertial  and  elastic 
structural  forces,  aerodynamic  forces  induced  by  the  static  or  dynamic  deformation  of  the 
structure,  and  external  disturbance  forces.  The  equation  of  motion  of  the  aeroelastic  system  in 
terms  of  discrete  system  can  be  derived  based  on  the  equilibrium  condition  of  these  forces,  i.e.. 

Mx(t)  +  Kx(t)  =  F(0  (2.1) 

where  M  and  K  are  the  mass  and  stiffness  matrices  generated  by  the  structural  finite  element 
method.  x(0  is  the  structural  deformation. 

The  structural  damping  matrix  is  excluded  in  Eq  2.1  for  simplicity,  but  it  can  be  easily  included. 
In  Eq  2.1,  the  terms  Mx(0  and  Kx(r)  are  the  inertial  and  elastic  structural  forces, 
respectively,  whereas  F(/)  represents  the  aerodynamic  forces  applied  on  the  structure.  F(0  can 
be  generally  split  into  two  parts;  the  aerodynamic  forces  induced  by  the  structural  deformation 
Fg  (x)  and  the  external  forces  F^ (t) ,  i.e.: 

F(0  =  F,(x)  +  F.(0  (2.2) 

The  external  forces  Fj(0  are  usually  provided.  Typical  example  of  Fj(r)  is  the  continuous 
atmospheric  turbulence  or  impulsive-type  gusts.  The  generation  of  Fa(x)  normally  relies  on  the 
theoretical  prediction  that  requires  the  unsteady  aerodynamic  computations.  Since  F^Cx) 

depends  on  the  structural  deformation  x(0,  the  relationship  can  be  interpreted  as  an  aerodynamic 
feedback.  Fig  2.1  presents  a  functional  diagram  that  illustrates  the  aeroelastic  interaction  of 
these  structural  and  aerodynamic  forces. 
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Fig  2.1  Aeroelastic  Functional  Diagram 

Without  the  aerodynamic  feedback,  Fig  2.1  reduces  to  an  open-loop  forced-structural  vibration 
system  whose  response  amplitude  is  usually  finite.  With  the  inclusion  of  F,(x),  Fig  2.1 

represents  a  closed-loop  dynamic  response  problem  which  can  be  expressed  by  the  following 
equation: 

Mx(r)  +  Kx(0-F,(x)  =  F.(0  (2.3) 

Eq  2.3  is  obtained  by  combining  Eq  2.1  and  Eq  2.2.  The  left  hand  side  of  Eq  2.3  is  in  fact  a 
closed-loop  aeroelastic  system  which  can  be  self-excited  in  nature.  This  gives  rise  a  stability 
problem  of  the  closed-loop  aeroelastic  system  known  as  flutter.  Flutter  analysis  usually  involves 
the  search  of  the  structural  stability  boundary  of  an  aircraft  structure  in  terms  of  its  flight  speed 
and  altitude  or  the  corresponding  dynamic  pressure.  If  F^  (x)  is  a  nonlinear  function  with  respect 
to  x(0,  the  flutter  analysis  must  be  performed  by  a  time-marching  procedure  solving  the 
following  equation: 

Mx(0  +  Kx(0  -  F.(x)  =  0  (2.4) 

with  initial  condition  of  x(0)  and  x(0)  being  specified  at  /  =  0. 

The  stability  boundary  of  the  aeroelastic  system  is  then  determined  by  examining  the  decay  or 
growth  of  the  structural  response  x(0  with  respect  to  the  flight  speed.  This  time-marching 
computational  procedure  is  rather  costly  since  it  generally  requires  a  nonlinear  time-domain 
unsteady  aerodynamic  method  known  as  the  Computational  Fluid  Dynamics  (CFD)  method. 
However,  the  conventional  industrial  practice  of  flutter  analysis  is  to  recast  Eq  2.4  into  a  set  of 
linear  systems  and  to  determine  the  flutter  boxmdaiy  by  solving  the  complex  eigenvalues  of  the 
linear  systems.  Such  a  procedure  involves  the  assumption  of  amplitude  linearization.  The 
amplitude  linearization  states  that  the  aerodynamic  response  varies  linearly  with  respect  to  the 
amplitude  of  the  structure  deformation  in  a  given  aircraft  motion  if  the  amplitude  is  sufficiently 
small  at  all  times. 

Since  flutter  analysis  is  a  dynamic  aeroelastic  stability  problem,  the  required  amplitude  for 
determining  such  a  stability  boundary  can  be  assumed  to  be  mathematically  infinitesimal.  This 
follows  that  the  amplitude  linearization  assumption  could  recast  Eq  2.4  into  an  eigenvalue 
problem  for  flutter  analysis.  In  this  case,  the  aerodynamic  system  can  be  approximated  by  a 
linear  system  for  which  an  aerodynamic  transfer  function  can  ^  defined.  This  transfer  function 
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relates  the  aerodynamic  feedback  F,(x)  to  the  structural  deformation  x(t)  by  means  of  the 
following  convolution  integral,  i.e.; 

F,(x)  = 

where: 

qcoB.  represents  the  aerodynamic  transfer  function 

^00  is  the  dynamic  pressure. 

L  is  the  reference  length  and  is  generally  defined  as: 

L  -  —  where  c  is  the  reference  chord 
2 

and: 

V  is  the  velocity  of  undisturbed  flow. 

The  Laplace  domain  counterpart  of  Eq  2.5  is  simply: 

Fa  (x(^))  =  ^=0  H  (  ^ )  x(^)  (2.6) 

where: 

H  is  the  Laplace  domain  counterpart  of  H 

With  Eq  2.6  at  hand,  Eq  2.4  can  be  readily  transformed  into  the  Laplace  domain  and  results  in  an 
eigenv^ue  problem  in  terms  of  s.  This  reads: 

[^5"M  +  K-^„H(^)jx(5)  =  0  (2.7) 

Since  the  finite  element  model  of  aircraft  structure  normally  contains  a  large  amount  of  degrees 
of  freedom,  the  size  of  the  mass  and  stiffness  matrices  are  usually  very  large.  Hence,  solving  the 
eigenvalue  problem  of  Eq  2.7  directly  would  be  computationally  costly.  To  circumvent  this 
problem,  one  introduces  the  “modal  approach”  which  can  be  expressed  as: 

X  =  <E>  q  (2.8) 

where  <(>  is  the  modal  matrix  whose  columns  contain  the  lower  order  natural  modes.  Normally, 
no  more  than  ten  numbers  of  the  lowest  natural  modes  are  sufficient  for  the  flutter  analysis  of  a 
wing  structure.  For  the  whole  aircraft  structure,  fifty  natural  modes  are  usually  sufficient,  q  is  the 
so-called  generalized  coordinates  which  are  the  eigenvectors  to  be  determined. 
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The  rationale  of  the  modal  approach  is  based  on  the  premises  that  the  critical  flutter  modes  are 
usually  due  to  the  coupling  of  the  lower  order  structural  modes.  Thus,  the  structural  deformation 
of  the  flutter  mode  can  be  sufficiently  represented  by  the  superposition  of  lower  order  modes. 

Substituting  Eq  2.8  into  Eq  2.7  and  pre-multiplying  Eq  2.7  withO^  yield  the  dynamic  (or  the 
flutter)  equation: 


where: 


and: 


q  =  0 

(2.9) 

M  =  M  O 

is  the  generalized  mass  matrix 

(2.10) 

K  =  K  $ 

is  the  generalized  stiffness  matrix 

(2.11) 

is  the  generalized  aerodynamic  forces  matrix. 

(2.12) 

Since  the  size  of  the  matrices  in  Eq  2.9  is  in  the  order  of  number  of  modes,  the  modal  approach 
provide  a  reduced  size  of  the  eigenvalue  problem  for  the  flutter  analysis  that  is  expressed  in  the 
generalized  coordinates.  Solving  such  a  reduced  size  eigenvalue  problem  is  computationally 
much  more  efficient  than  that  of  Eq  2.7.  Eq  2.9  is  generally  referred  as  the  classical  flutter 
matrix  equation. 

The  above  discussion  shows  that  reducing  the  time-domain,  generally  non-linear  flutter  equation 
(Eq  2.4)  to  the  classical  flutter  matrix  equation  (Eq  2.9)  lies  in  the  availability  of  the 
aerodynamic  transfer  function.  However,  the  generation  of  aerodynamic  transfer  functions  in  the 
Laplace  domain  by  solving  unsteady  aerodynamics  can  be  a  veiy  complicated  procedure.  For 
this  reason,  imsteady  aerodynamic  methods  are  often  formulated  in  the  frequency  domain  by 
assuming  simple  harmonic  motion.  The  frequency-domain  aerodynamic  transfer  function  in 
matrix  form  is  called  the  Aerodynamic  Influence  Coefficient  (AIC)  matrix.  In  fact,  the  major 
functionality  of  ZAERO  is  to  generate  such  AIC  matrices  for  the  aircraft  configuration.  This  will 
be  briefly  discussed  in  the  following  section. 

2.2  Unified  MC  of  ZAERO 

Four  unsteady  aerodynamic  methods  are  incorporated  in  ZAERO,  namely  ZONA6,  ZONA7, 
ZTAIC  and  ZONA7U  that  jointly  generate  the  AIC  matrices  covering  the  complete  domain  of 
Mach  number  range.  Theses  ZAERO  generated  AIC  matrices  are  called  the  mified  AIC  for  their 
unified  feature  in  Mach  number  range. 

ZONA6  and  ZONA7  solve  the  integral  equations  due  to  the  respective  unsteady  linearized  small- 
disturbance  subsonic  and  supersonic  equations  for  general  aircraft  configurations.  The  integral 
equation  is  formulated  in  the  reduced  frequency  domain  which  is  in  the  context  of  the  simple 
harmonic  motion.  The  reduced  frequency,  denoted  as  k,  is  a  fundamental  unsteady  aerodynamic 
parameter  defined  as: 
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(2.13) 


k  = 


coL 

~V~ 


where: 

Of  is  the  harmonic  oscillatoiy  frequency. 

Panel  method  is  adopted  by  ZONA6  and  ZONA7  to  solve  the  integral  equations.  The  panel 
method  requires  the  aircraft  configuration  to  be  discretized  into  many  small  panels.  Each  panel 
is  defined  as  “aerodynamic  box”.  Fig  2.2  shows  a  typical  panel  model  of  a  wing-body 
configuration. 


Fig  2.2  Typical  panel  model  of  a  wing-body  configuration 

Each  aerodynamic  box  contains  a  control  point  where  the  boundary  condition  is  imposed. 
According  to  the  panel  model,  the  integral  equation  is  then  approximated  by  the  summarization 
of  elementary  integrals  associated  with  each  aerodynamic  boxes.  The  assembly  of  the 
elementary  integral  solutions  gives  a  matrix  whose  coefficients  represent  the  aerodynamic 
influence  of  aerodynamic  boxes  to  control  points.  This  matrix  is  called  the  Aerodynamic 
Influence  Coefficient  (AlC)  matrix  which  relates  the  structural  deformation  to  the  aerodynamic 
forces  by: 

(2.14) 

where: 

h  is  the  structural  deformation  defined  at  the  aerodynamic  boxes 

Fh  is  the  resultant  aerodynamic  forces  at  the  aerodynamic  boxes  due  to  h. 
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Comparing  Eq  2.4  with  Eq  2.6,  one  can  see  that  [  AlC(ik)  ]  is  indeed  the  aerodynamic 
transfer  function,  except  that: 

•  AIC  matrix  is  computed  in  the  reduced-frequency  domain  (^-domain)  rather  than  the 
Laplace  domain. 

•  AIC  matrix  is  computed  based  on  the  panel  model  which  is,  in  general,  considerably 
different  from  its  respective  structural  finite  element  model.  This  gives  rise  the 
problem  of  displacement  and  force  transferal  between  these  two  models.  These 
issues  can  be  resolved  by  a  spline  matrix  that  relates  or  interpolates  the 
displacements  at  structural  finite  element  grid  points  to  those  at  aerodynamic  panel 
model.  The  derivation  of  spline  matrix  will  be  discussed  later.  The  detailed 
theoretical  background  of  ZONA6  and  ZONA7  is  discussed  in  Section  3. 

The  ZTAIC  method  computes  the  unsteady  aerodynamics  for  transonic  flows.  It  requires  the 
steady  pressure  distribution  input  from  the  users.  ZTAIC  utilizes  an  inverse  design  method  for 
airfoil  design  by  solving  the  transonic  small-disturbance  equation.  These  airfoil  sections  produce 
the  steady  pressure  distribution  that  matches  input  steady  pressure.  Unsteady  pressure 
coefficients  on  the  airfoil  sections  are  then  computed  by  solving  the  unsteady  transonic  small 
disturbance  equation.  These  unsteady  pressure  coefficients  are  used  to  correct  the  linear 
unsteady  pressure  computed  by  ZONA6  for  the  transonic  shock  effects.  At  the  end,  ZTAIC 
employs  a  modal  AIC  approach  to  reconstruct  the  AIC  matrix  computed  by  ZONA6,  leading  to  a 
transonic  AIC  matrix.  For  detailed  theoretical  description,  see  Section  4  for  the  formulation  of 
the  ZTAIC  method. 

ZONA7U  is  a  unified  supersonic/hypersonic  unsteady  aerodynamic  method.  This  method 
combines  the  supersonic  lifting  surface  aerodynamics  of  ZONA7  with  a  nonlinear  correction 
matrix  based  on  Donov  &  LinnelTs  imiformly-valid,  hi^-order  Hypersonic/Supersonic  Scheme. 
This  correction  matrix  takes  the  flow  nonlinearity  as  well  as  the  flow  rotationality  due  to 
oscillatory  shock  waves  into  an  account,  which  covers  both  the  Mach  wave  and  Newtonian 
limits.  The  only  additional  input  for  ZONA7U  is  the  sectional  airfoil  thickness  distribution. 
ZONA7U  generates  the  AIC  matrix  that  is  in  the  same  form  of  that  of  ZONA7. 

The  detailed  theoretical  formulations  of  ZONA7U  is  discussed  in  Section  5. 

2.3  Functionality  of  the  Spline  Matrix 

The  problem  of  data  transferal  between  the  panel  model  and  the  structural  finite  element  model 
usually  amounts  to  the  displacement  transferal  from  the  structural  grid  points  to  the  aerodynamic 
control  points  of  the  panel  model  and  that  of  the  forces  from  the  aerodynamic  control  points  to 
structural  grid  points.  ZAERO  provides  a  spline  module  which  generates  a  spline  matrix  G  such 
that: 


h  =  Gx  (2.15) 

Note  that  h  is  defined  by  the  aerodynamic  control  points  whereas  x  by  the  structural  finite 
element  grid  points.  Four  methods  are  incorporated  in  the  spline  module,  namely  the  beam 
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spline  method,  infinite  plate  spline  method,  thin-plate  spline  method  and  rigid-body  attachment 
method.  These  four  spline  methods  jointly  construct  the  spline  matrix  G  for  the  displacement 
transferal.  The  detail  formulations  of  these  four  spline  methods  are  discussed  in  Section  6. 

Once  the  spline  matrix  G  is  generated,  the  force  transferal  from  the  aerodynamic  control  points 
to  structural  grid  points  can  be  performed  accordingly: 

F,  =  G^  (2.16) 

Eq  2. 16  is  derived  based  on  the  principle  of  virtual  work.  Since  the  forces  at  aerodynamic  boxes 
Fh  and  their  structurally  equivalent  values  Fa  must  do  the  same  virtual  work  in  their  respective 
displacements,  we  have: 

Fj,  =  F^  (2.17) 


where: 

6h  and  5x  are  virtual  displacements. 

Substituting  Eq  2. 1 5  into  the  left-hand  side  of  Eq  2. 17  and  upon  rearranging,  yields: 

(f,  -  G^  Ft  )  =  0  (2.18) 

Because  of  the  arbitrariness  of  the  virtual  deflection,  i.e.  6x  ^  0,  the  terms  in  the  bracket  of  Eq 
2.18  must  vanish,  which  leads  to  Eq  2.16.  With  Eq  2.15  andEq  2.16  at  hand,  the  aerodynamic 
feedback  acting  at  the  structural  grid  points  can  be  obtained.  Combining  Eq  2.14  and  Eq  2.15 
and  substituting  the  resultant  equation  into  Eq  2. 16  yield: 

F,  =  q^G'^  [AIC{ik)]G  x  (2.19) 

Again,  note  the  similarity  between  Eq  2.19  and  Eq  2.6.  The  generalized  aerodynamic  forces 
matrix  in  the  it-domain  can  be  obtained  and  computed  by  applying  the  modal  approach: 

Q{ik)  =  G”"  [  AlCiik)  ]  G  O  (2.20) 

The  it-domain  Q(/A)  results  from  the  reduction  to  the  simple  harmonic  motion  (frequency 
domain)  from  the  transient  motion  (Laplace  domain)  introduced  previously  in  the  unsteady 
aerodynamic  formulation.  The  impact  of  ^-domain  Q(/^)  on  the  solution  technique  required  to 
solve  the  classical  flutter  matrix  equation  will  be  discussed  in  the  following  section: 

2.4  Impact  of  Q(iA:)  on  Flutter  Solution  Technique 

The  reduction  to  the  harmonic  motion  introduced  in  the  previous  unsteady  aerodynamic 
formulation  indicates  that  the  ifc-domain  Q(/^)  is  valid  only  for  steady  state  response  of  the 
structure.  This  also  implies  that  the  flutter  solution  is  valid  only  at  the  flutter  boundary  where 
damping  of  the  aeroelastic  system  is  zero.  Since  the  generalized  aerodynamic  forces  matrix  is 
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available  only  in  the  frequency  domain,  the  frequency-domain  counterpart  of  Eq  2.9  can  be 

s  L 

obtained  by  replacing  Q(— )  byQ(/A:)  andj  by  m.  This  gives: 

+  Q(ik)]  q  =  0  (2.21) 

As  a  stability  measure,  an  artificial  structural  damping  is  added  to  Eq  2.21: 

[-u>^  M  +  (1  +  /gJK  -  Q(/^)]q  =  0  (2.22) 


where: 

gs  is  the  added  artificial  structural  damping. 

Eq  2.22  is  the  so-called  AT-method  flutter  equation.  It  is  mathematically  consistent  with  the 
assumption  of  simple  harmonic  motion  introduced  in  the  unsteady  aerodynamic  formulation. 
But  its  predicted  damping  is  only  an  artifice  used  to  seek  out  the  flutter  point  and  cannot  be 
interpreted  as  having  a  physical  significance  as  a  measure  of  decay  rate  of  the  aeroelastic 
response. 

While  the  AT-method  flutter  equation  is  capable  providing  the  prediction  of  the  flutter  boundary, 
e.g.  flutter  speed  at  zero  damping,  it  is  often  desired  to  have  a  reliable  damping  prediction  to 
detect  the  aeroelastic  characteristics  at  sub-critical  flight  speeds.  The  predicted  damping  values  at 
sub-critical  flight  speeds  can  serve  as  a  guideline  for  conducting  wind  tunnel  or  flight  flutter 
tests.  It  is  for  this  reason  that  the  P-K  method  is  widely  adopted  by  aeroelasticians  as  the 
primary  tool  for  finding  flutter  solutions.  In  Eq  2.9,  the  P-K  method  replaces  the  Laplace 

sL 

domain  generalized  aerodynamic  forces  matrix  Q(— )  by  Q(/^);  and  it  further  defines  a  non- 
dimensional  Laplace  parameter p  such  that: 

p  =  ^  =  yc  +  ik)  (2.23) 


where: 

y  is  the  transient  decay  rate  coefficient, 

In  this  way,  Eq  2.9  becomes: 


q  =  0 


(2.24) 


For  a  given  Q(/ifc),  the  P-K  method  solves  the  complex  eigenvalues  of  Eq  2.24  in  terms  of  p.  Eq 
2.24  is  solved  at  several  given  values  of  V  and  q^  ,  for  complex  eigenvalues  p  associated  with 
modes  of  interest.  This  is  accomplished  by  an  iterative  procedure  that  matches  the  reduced 
frequency  k  to  the  imaginary  part  of  p  for  every  structural  mode. 
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Strictly  speaking,  the  P-K  method  is  mathematically  inconsistent  because  Q(/A:)  is  obtained  by 
the  simplification  of  simple  harmonic  motion  (or  by  applying  the  Fourier  Transform)  but  other 
terms  in  Eq  2.24  are  associated  with  the  transient  motion  (or  by  applying  the  Laplace 
Transform).  However,  it  is  generally  believed  that  the  P~K  method  is  a  good  approximation 
method  in  finding  the  rate-of-decay  type  of  solution.  The  rationale  for  the  P-K  method  is  that 
for  simple  harmonic  motion  with  slowly  increasing  or  decreasing  amplitude,  the  generalized 
aerodynamic  force  based  on  constant  amplitude  are  good  approximations. 

The  detailed  formulations  and  the  solution  technique  of  the  AT-method  and  the  P-K  method  will 
be  discussed  in  Section  7.  In  Section  7,  a  newly  developed  flutter  solution  method  developed  by 
ZONA  Technology  called  the  g-method  will  also  be  presented.  The  g-method  includes  a  first 
order  damping  term  in  the  flutter  equation  that  is  rigorously  derived  from  the  Laplace-domain 
aerodynamics.  The  g-method  generalizes  the  .K^-method  and  the  P-K  method  and  could  provide 
unlimited  roots  of  the  flutter  equation.  The  extra  roots  obtained  by  the  g-method  are  associated 
with  the  aerodynamic  lag  root  which  could  not  be  otherwise  obtained  by  the  two  methods.  The 
superiority  and  the  physical  significance  of  the  g-method  will  also  be  discussed  in  Section  7. 
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3.0  ZONA6/ZONA7:  UNSTEADY  SUBSONIC  /  SUPERSONIC 
AERODYNAMICS  FOR  WING-BODY  AIRCRAFT 
CONFIGURATIONS  WITH  EXTERNAL  STORES 


ZONA6  and  ZONA7  solve  the  respective  unsteady  three-dimensional  linearized  small- 
disturbance  potential  equations  of  subsonic  and  supersonic  aerodynamics.  Their  unique  feature 
lies  in  the  capability  of  modeling  realistic  configurations  such  as  an  aircrafl/wing-body 
combination  including  external  stores  or  nacelles.  In  this  section,  we  will  discuss: 

•  Background  of  ZONA6  and  ZONA  7. 

•  Derivations  of  unsteady  boundary  condition  and  pressure  coefficient. 

•  Derivations  of  the  unsteady  linearized  small-disturbance  equations  and  their  corresponding 
integral  equations. 

•  Formulation  of  panel  method  employed  by  ZONA6  and  ZONA7  for  solving  the  integral 
equations. 

•  Generation  of  their  Aerodynamic  Influence  Coefficient  (AIC)  matrices. 

•  Treatments  of  body-wake  effect  in  subsonic  flow. 

•  Treatments  of  inlet  boxes  and  superinclined  aerodynamic  boxes  in  supersonic  flow. 

3.1  Backgrounds  of  ZONA6  and  ZONA7 

Since  1985,  ZONA  has  been  devoting  its  R&D  effort  to  the  development  of  unsteady 
aerodynamic  methods  for  aeroelastic  applications.  The  first  ZONA  software  product  for 
supersonic  lifting  surface  unsteady  aerodynamics  is  the  ZONA51  code.  ZONA51  employs  the 
acceleration-potential  approach  for  thin-wing  type  of  lifting  surfaces  (Ref  3.1).  This 
acceleration-potential  approach  is  the  outgrowth  from  the  Harmonic  Gradient  Method  (HGM) 
developed  by  Chen  and  Liu  in  1985  (Ref  3.2).  Today,  ZONA51  is  the  industrial  standard 
method  for  supersonic  lifting  surface  unsteady  aerodynamics  in  MSC/NASTRAN  -  Aero  Option 
II  (Ref  3.3). 

In  Figs  3.1  and  3.2,  computed  results  of  Re(ACp)  and  Im(ACp)  of  a  flat  plate  undergoing 

plunging  and  pitching  motion  about  the  leading  edge  are  presented.  It  can  be  seen  that  the  results 
computed  by  ZONA5 1  are  in  better  agreement  with  Jordan’s  exact  solution  (Ref  3.4)  than  that  of 
HGM.  The  close  agreement  with  Jordan’s  result  near  the  trailing  edge  is  attributed  to  an  exact 
treatment  of  linearized  pressure  in  the  acceleration-potential  formulation  of  ZONA51. 

At  the  National  Aeronautical  Establishment  (NAER)  Canada,  Lee  (Ref  3.5)  has  performed  a 
comparative  study  on  computed  pressures  for  an  oscillating  leading-edge  flap  on  the  F-18  wing 
at  58.8%  span  and  at  M=  1.1  and  reduced  frequency  k  =  4.0  (see  Fig  3.3).  Two  computer  codes 
are  used:  the  ZONA51  code  and  the  SPBP  code  (Ref  3.6). 
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Fig  3.1  Comparison  of  Computed  Unsteady 
Pressures  for  a  Plunging  M  - 1.25  and  ft  =  2.0 


In  Fig  3.4,  computed  results  show  that  the 
expected  pressure  jump  at  the  hinge  line  is 
well  predicted  by  ZONA51  using  126 
panels.  The  results  obtained  by  the  SPP 
code,  however,  depend  on  the  number  of 
quadrature  points  chosen  where  no 
pressure  jump  across  the  hinge  line  is  seen 
until  more  quadrature  points  are  used. 

Johnson  et  al.  (Ref  3.7)  performed  as 
supersonic  unsteady  aerodynamic 
computation  for  the  Viggen  Idealization 
(Ref  3.8)  using  ZONA5 1  option  in 
MSC/NASTRAN.  The  canard  is  slightly 
above  the  vsing  plane  by  a  distance  of  0. 1 
length  units.  The  generalized  force  Q12 
shown  is  the  lift  coefficient  of  the  wing 
due  to  a  unit  rotation  of  the  canard  about 
its  midchord. 


M  *  1.25.  k  ■  2.0 


Fig  3.2  Comparison  of  Computed  Unsteady 
Pressures  for  a  Pitching  Flat  Plate  About  the 
Leading  Edge  at  Af  =  1.25  and  k  =  2.0 


j  IQM 


Fig  3.3  Aerodynamic  Modeling  of  F-18  Wing 


The  lift  coefficient  is  obtained  by  dividing  die  generalized  force  by  the  dynamic  pressure  and  the 
square  of  the  wing  span  (1.20  units).  As  shown  in  Ref  3.6,  the  wing  lift  is  calculated  by 
considering  two  rigid-body  modes  of  motion,  the  first  being  the  wing  alone  undergoing  unit 
plunging,  and  the  second  being  unit  canard  rotation.  The  reduced  frequency  is  based  on  a 
reference  semichord  of  1.00  unit.  The  results  are  shown  in  Fig  3.5  and  almost  coincide  with  the 
NLR  SPNLRI-CP  (Ref  3.9)  results  up  to  the  high  reduced  frequency  ofk  =  5.0. 
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To  verify  convergence,  two  aerodynamic  element  (“box”)  idealizations  were  considered.  In  the 
first,  the  canard  was  divided  into  8  equal  width  strips  with  8  equal  chordwise  divisions,  and  the 
wing  was  divided  into  12  equal  width  strips  and  12  equal  chordwise  divisions,  giving  a  total  of 
208  boxes.  In  the  second  idealization,  the  equal  divisions  were  10  x  10  on  the  canard  and  20  x 
20  on  the  wing  for  a  total  of  500  boxes.  The  two  sets  of  results  agreed  within  plotting  accuracy 
and  are  not  distinguished  in  Fig  3.5. 

Stark’s  Characteristic  Box  Method  (CHB;  Ref  3.10)  results  are  also  shown  in  Fig  3.5.  It  can  be 
seen  that  the  CHB  results  agree  with  the  NLR  SPNLRI-CP  and  ZONA51  results  only  at  low 
reduced  fi-equencies  (^  <  1).  At  higher  reduced  frequencies,  the  CHB  results  give  considerable 
discrepancy  compared  to  these  of  the  other  two  methods. 


Fig  3.4  Effects  on  Panels  on  Computed  Pressures 
for  an  F-IS  Wing  with  an  Oscillating  Leading-Edge 
Flap  at  58.8%  span, M=  1.1  and  k  =  4.0 


Fig  3.5  Real  and  Imaginary  Part  of  the  Lift  on  the 
Main  Wing  Due  to  Pitch  of  the  Canard,  M  =  1.054 
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In  order  to  demonstrate  the  capability  of 
ZONA51  for  supersonic  flutter  analysis,  a  70® 
swept  delta  wing  (model  1  A)  in  this  series  is 
selected  for  comparison  with  the  existing  flutter 
results.  The  planform  is  subdivided  into  100 
panels  as  shown  in  Fig  3.6.  Four  modes  are 
used  in  the  present  flutter  analysis. 


Fig  3.6  Paneling  Scheme  for  a  70®  Delta  Wing 


A  set  of  flutter  points  were  obtained  for  seven  Mach  numbers  (A/=  1.01,  1.19,  1.30,  1.64,  2.0, 
2.25  and  3.0)  using  ZONA51.  Throughout  the  supersonic  Mach  number  range  considered, 
ZONA51  appears  to  yield  the  best  correlation  with  experimental  data  among  all  supersonic 
methods  shown  here;  these  include  FAST  (Ref  3. 1 1),  ACUNN  (Ref  3. 12),  Piston  Theory  (second 
order)  (Ref  3.13)  and  CAP-TSD  (Ref  3.14)  (Figs  3.7,  3.8).  Overall,  ZONA51  predicts  slightly 
lower  values  in  the  flutter  frequency  ratio  Of  /Oj  (®2  the  natural  frequency  of  mode  2)  than 
the  obtained  experimental  data. 


O  CAP-TSD 
o  PMonTtMory 
A  ZONA51C 
A  FAST 


M 


Fig  3.7  Computed  and  Measured  Flutter  Speeds 
vs.  Mach  Number 


Fig  3.8  Comparison  of  Ratios  of  Flutter  Frequency 
vs.  Mach  Number 


ZONA7  (Ref  3.15)  generalizes  ZONA51  for  the  wing-body  configuration.  Its  lifting  surface 
method  is  identical  to  ZONA51.  But  its  body  aerodynamic  capability  enables  ZONA7  to  model 
realistic  aircraft  configuration  including  the  external  stores.  To  verify  ZONA7  on  body-alone 
configuration,  a  10%  thick  ogive  body  performing  oscillation  in  a  first  bending  model  is  selected 
(Fig  3.9).  Complete  agreement  is  found  between  the  HPP  (Ref  3.16)  results  and  the  present 
results  in  this  hi^  reduced-frequency  case  {k=  1.0). 
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Fig  3.9  In-Phase  and  Out-Of-Phase  Pressure  Coefficients  for  a  Bolic-Ogive  in  First  Bending  Mode 

at  Af«o  =  2.5  and  Reduced  Frequency  k  =  1.0 

For  the  case  of  wing-body  configuration  in  steady  flow,  longitudinal  loadings  {C^dld^^ )  over  a 
10%  thick  body  with  and  without  a  tapered  wing  (AR  =  4.0  and  taper  ratio  =  0.6)  are  presented  in 
Fig  3.10. 


PERCENT  FUSELAGE  LENGTH 


Fig  3.10  Static  Loadings  About  the  Fuselage-Wing  Junctions 
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No  interference  effect  is  noted  in  the  first  40%  of  the  body  length,  as  expected.  A  bump-like 
loading  along  40-85%  body  length  is  observed  as  a  result  of  the  presence  of  the  tapered  wing.  In 
addition,  lifting  pressure  distribution  along  the  wing  chord  is  also  plotted  at  the  20% 
semispanwise  location.  Good  agreement  is  found  between  the  computed  and  measured  results 
(Ref  3.17).  Thus,  the  steady  aerodynamic  option  of  the  ZONA7  method  is  validated  by  means  of 
this  present  interference  example. 


For  wing-store  configuration,  an  NLR  wind- 
tunnel  test  configuration  constructed  with  an 
F-5  wing  plus  an  underwing  store  is  modeled 
by  112  body  panels  representing  the  missile 
body,  72  wing  panels  for  the  launcher  and  24 
panels  for  the  four  aft  fins  (Fig  3.11).  The 
complete  configuration  is  in  pitching 
oscillation  about  50%  root  chord  at  a  circular 
frequency  f  —  10  Hz,  and  at  two  Mach 
numbers,  M=  1.1  and  1.35. 


Fig  3.11  Paneling  Model  for  the  Underwing  Store 
Configuration:  Northrop  F-5  Wing  Plus  Underwing 
Pylon,  Launcher,  Missile  Body  with  Four  Aft  Fins 


Figs  3.12-3.14  present  the  comparisons  of  the  NLR  measured  data  and  the  present  computed 
results  are  presented. 


Fig  3.12  presents  the  unsteady 
normal  forces  and  pitching 
moments  on  the  imderwing  store 
system  of  two  different 
combinations:  1)  the  pylon  plus 
laimcher  {P  +  L),  and  2)  in 
addition  to  1),  a  missile  body 
with  four  aft  fins  (P  +  L  +  MB  + 
AW).  The  in-phase  (real)  part  of 
the  computed  normal  forces  and 
pitching  moments  for  both  cases 
correlated  well  with  the 
measured  data,  while  all  out-of¬ 
phase  (imaginary)  parts  remain 
insignificantly  small. 


TEST  DATA  (NLR)-  F*20Hx 
CONFIGURATION 
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Fig  3*12  Unsteady  Normal  Force  and  Pitching  Moment  for  the 
Underwing  Store  Configuration  With  and  Without  the  Missile  Body 
Oscillating  about  50%  Root  Chord  at  Mao  =  1.1  and  1.35  and  Reduced 
Frequency  k  ==  0.1 
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Fig  3.13  Unsteady  Side  Force  and  Yawing  Moment 
for  the  Complete  Underwing  Store  Configuration 
With  and  Without  the  Missile  Body  and  Launcher 
Oscillating  about  50%  Root  Chord  atM^,  =  1.1  and 
1.35  and  Reduced  Frequency  k  =  0.1 


Fig  3.14  Unsteady  Spanwise  Normal  Force  and 
Pitching  Moment  for  the  Clean  F-5  Wing  and  the 
Underwing  Store  Configuration  at  =  1.35  and 
Reduced  Frequency  k  =  0.1 


Fig  3.13  presents  the  side  forces  and  yawing  moments  on  the  underwing  store  systems  in  three 
different  combinations:  1)  pylon  alone  (P);  2)  pylon  and  launcher  {P  +  P);  and  3)  in  addition  to 
2),  a  missile  body  with  four  aft  fins  (P  +  L  +  AfB  +  AW).  Both  measured  data  and  computed 
results  show  increases  in  the  in-phase  normal  forces  (positive  for  outboard  direction)  with  the 
addition  of  the  system  from  1)  to  3),  whereas  a  decreasing  trend  is  observed  for  the  out-of-phase 
parts.  The  unsteady  yawing  moments  on  the  system  are  relatively  small.  The  computed  results 
also  show  that  the  added  missile  with  fins  to  the  system  contributes  most  to  the  in-phase  moment 
(positive  body  apex  pointing  inboard). 

In  Fig  3.14,  the  integrated  spanwise  unsteady  normal  forces  and  pitching  moments  along  the  F-5 
wing  under  the  influence  of  the  complete  underwing  store  system  are  plotted  against  those  of  the 
clean-wing  case  according  to  the  computed  results  and  the  test  data.  It  seen  that  the  computed 
forces  and  moments  predict  the  same  trend  as  the  measured  data  showing  a  finite  discontinuity 
across  the  pylon  location.  The  computed  results  tend  to  overestimate  the  in-phase  forces  and 
moments  and  underestimate  the  out-of-phase  forces  and  moments  in  comparison  with  the 
measured  data. 


However,  these  discrepancies  may  be  caused  by  the  uncertainties  in  the  measured  unsteady  data 
as  mentioned  in  Ref  3.18.  Meanwhile,  NLR  also  provided  the  measured  quasisteady  data,  which 
is  supposedly  more  reliable  for  the  in-phase  forces.  Better  agreement  in  trend  between  the 
quasisteady  data  and  the  computed  results  is  found  for  this  case. 
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Z0NA6  (Ref  3.19)  is  the  subsonic  counterpart  of  ZONA7  except  that  includes  the  important 
body-wake  effects  for  fuselage  and  stores.  For  wing-body  configuration  in  steady  flow,  the  static 
loading  of  the  same  fuselage-wing  configuration  depicted  in  Fig  3.10  is  now  computed  by 
ZONA6  at  A/  =  0.6,  a  =  4°.  Fig  3.15  shows  an  excellent  agreement  found  between  ZONA6 
results  and  the  measured  data. 


Ptretnt  Fustl»gt  Ltngth 


F^3.15  Static  Loading  of  NACA  Wing-Body  Configuration;  jlf=0.6anda  =  4® 

Figs  3.16  -  3.19  present  a  series  of  computed  results  by  ZONA6  in  comparison  with  NLR 
measured  and  analysis  data  for  an  NLR  modeled  wing-tip-store  configuration  (Ref  3.20). 


Based  on  the  configuration  and  the  paneling 
scheme  shown  in  Fig  3.16,  the  computed 
steady  pressure  coefficients  along  the  tip- 
tank  at  two  azimuthal  angles,  1)  0  =  292.5° 
and  2)  6  =  67.5°,  are  presented  in  Fig  3.17. 
The  present  results  appear  to  be  in  good 
agreement  with  each  other  and  with  NLR’s 
computed  and  measured  data. 


Fig  3.16  NLR  Wing-Tip-Tank  Configuration  Showing 
Paneling  Scheme 


Fig  3.18  presents  the  unsteady  pressures  along  tip-tank  at  an  azimuthal  angle  0  =  202.5°.  It  is 
seen  that  the  present  results  are  in  good  agreement  with  the  NLR  measured  data.  The  present  in- 
phase  Cp  without  wake  appears  to  deteriorate  towards  the  tail  of  the  tip-tank,  whereas  that  with 
wake  appears  to  correlate  best  with  the  measured  data.  This  is  expected,  because  the  flow  is 
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known  to  separate  at  the  rear  in  this  case.  For  the  out-of-phase  Cp,  the  NLR  computed  results 
show  a  discrepancy  with  the  measured  data  starting  from  the  midbody. 


Fig  3.17  Steady  Pressure 
Distribution  Along  the  Tip-Tank 
of  NLR  Wing-Tip-Tank 
Configuration; 

M=  0.45  and  a  =  4“: 

a)  e  =  292.5“  and  5)0  =  67.5“ 


Fig  3.18  Unsteady  Pressure 
Distribution  Along  the  Tip- 
Tank  of  NLR  Wing-Tip-Tank 
Configuration; 

M=  0.45,  a = 0“,  A  =  0.305, 
X„=0.15cr,  and  0  =  202.5“ 


Fig  3.19  Unsteady  Normal  Load 
Distribution  Along  the  Tip-Tank 
of  NLR  Wing-Tip-Tank 
Configuration; 
M=0.45,a=0“,A  =  0305 
and  Xo  =  0.15  Cr, 


Fig  3.19  presents  the  unsteady  normal  load  Cz  along  the  tip-tank.  It  is  seen  that  the  present 
predicted  values  of  the  in-phase  and  out-of-phase  Cz  are  in  better  agreement  with  the  measured 
data  than  those  of  NLR,  particularly  for  the  cases  with  the  wake  model.  In  addition,  the  NLR’s 
out-of-phase  result  show  considerable  departure  from  the  measured  data. 

It  should  also  be  noted  that  ZONA6's  lifting  surface  method  adopts  a  higher  order  paneling 
scheme  than  the  Doublet  Lattice  Method  (DLM)  (Ref  3.21),  in  that  the  later  only  considers  a  line 
distribution  on  each  aerodynamic  box.  It  is  found  that  the  hi^  order  paneling  scheme  is  of 
importance  for  the  robustness  of  the  unsteady  lifting  surface  methods  (Ref  3.22).  In  order  to 
demonstrate  this,  two  wing  planforms  are  selected  for  comparison  of  the  present  results  and 
DLM  results.  These  planforms  include  a  rectangular  wing  of  aspect  ratio  20  at  A/=  0.0  and  a  70® 
delta  wing  at  M=  0.8.  The  main  objective  in  the  present  cases  studied  is  to  display  the  ease  of 
utilization  of  a  high-order  panel  code  such  as  ZONA6  as  opposed  to  that  of  DLM.  For  all  cases 
considered  (Figs  3.20-3.27),  a  commonly-practiced  paneling  scheme,  as  shown  in  Fig  3.22,  the 
spanwise  panel  cut  Ny  is  provided  by  equal  cut  along  the  span  whereas  the  chordwise  panel  cut 
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74  by  lines  that  connects  the  equally-divided  points  along  the  root  chord  and  tip  chord.  It  will  be 
shown  that  ZONA6  is  less  sensitive  to  the  paneling  scheme  than  DLM,  hence  a  more  robust 
code.  No  attempt  is  made  to  select  other  paneling  schemes  for  further  investigation  of  this  issue. 


(b)  10x40 


3.20  Generalized  Aerodynamic  Forces  Cl^  (=Qi2/S)  versus  Reduced  Frequency  k: 

AR  -  20,  Rectangular  Wing  (M  =  0.0,  Xo  =  0.5c); 
a)  Panel  Number  =  10  x  10,  b)  Panel  Number  =  10  x  40 

Rectaneular  Wine  of  AR  =  20 

Fig  3.20  presents  the  ZONA6  and  DLM  results  in  terms  of  generalized  aerodynamic  forces 
(GAF)  versus  reduced  frequencies  K  for  a  rectangular  wing  of  AR  =  20.0  pitching  about  50% 
chord.  Two  paneling  schemes  are  adopted:  i\4  x  =  10  x  10  (Fig  3.20a)  and  Nx.xNy—  10  x  40 
(Fig  3.20b).  It  can  be  seen  that  for  the  10  x  10  case  the  DLM  results  depart  from  the 
Theodorsen’s  exact  theory  (Ref  3.23)  exponentially  for  k  =  0.3,  while  the  ZONA6  results 
maintains  an  acceptable  range.  For  the  10  x  40  case,  the  error  of  ZONA6  results  are  further 
reduced  as  one  would  expect.  However,  large  discrepancy  still  exists  in  the  imaginary  part  of 
C,  between  DLM  and  Theodorsen’s  theory. 

According  to  the  MSC/NASTRAN  manual  (Ref  3.3),  the  chordwise  panel  number  is  linearly 
proportional  to  reduced  frequency  k.  For  obtaining  a  valid  result,  DLM  requires  that  k  =  0.08  i\4, 
for  at  least  equal  to  or  greater  than  four.  This  is  to  say  that  for  A4  =  10  the  valid  limit  for  k  is 
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around  1.0,  and  for  k  =  10,  has  to  be  over  100.  For  a  rectangular  wing,  the  spanwise  panel 
number  Ny  can  be  expressed  as: 


= 


(  AR 
l^AR 


\ 


where  AR  and  AAR  are  the  aspect  ratio  of  wing  and  panel,  respectively.  If  AAR  is  bounded  by 
3.0,  Ny  should  be  at  least  over  400.  It  is  clear  that  all  these  tedious  numerical  restrictions  are  a 
result  of  the  low  ordemess  of  DLM.  By  contrast,  ZONA6  removes  all  these  restrictions  and  can 
be  applied  readily  with  the  present  paneling  scheme.  A  further  study  of  these  methods  in  the 
effect  of  the  panel  number  on  and  results  is  shown  in  Fig  3.21.  It  is  seen  that  the  DLM 

result  is  very  sensitive  to  the  spanwise  panel  cut  if  Ay  is  below  20  and  it  converges  very  slowly  to 
Theodorsen’s  solution.  On  the  other  hand,  ZONA6  results  approach  Theodorsen’s  solution 
asymptotically  as  the  panel  number  increases. 


Fig  3.21  Effects  of  Panel  Numbers  on  GAP  and  ; 
AR  =  20;  Rectangular  Wing  (M  =  0.0,  k  =  10,  Xo  =  0.5c) 


Slender  Delta  Winss 


Fig  3.22  Typical  Paneling  Scheme 
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The  slender  delta  wing  family  belongs  to  another  planform  category  that,  with  the  present 
paneling  scheme  (see  Fig  3.22),  DLM  would  have  difficulty  in  producing  satisfactory  results. 
For  a  tapered  or  delta  wing  planform,  Ny  can  be  generalized  to: 


AAR,  = 


2K 


AR 


where  AAR,  is  the  panel  aspect  ratio  of  the  strip  and  c,-  is  the  corresponding  local  chord  length. 
Hence,  the  ratio  of  AAR,  on  two  strips  is  inversely  proportional  to  their  local  chord  lengths.  For 

a  70°  delta  wing  with  a  10  x  10  panel  cut  (Fig  3.22),  the  aspect  ratio  of  the  outboard  strip  at  the 
tip  section,  would  amount  to  17  times  that  of  an  inboard  strip  at  the  second  station.  This 
condition  would  probably  be  too  stringent  a  condition  for  DLM  to  be  applicable. 

The  effect  of  swept  angle  on  the  steady  and  for  a  slender  wing  is  presented  in  Fig  3.23 

at  M  =  0.8.  With  a  10  x  10  panel  cut  for  a  slender  delta  wing,  it  seen  that  ZONA6  result 
correctly  approaches  Miles’  asymptotic  slender-wing  limit  (Ref  3.24),  whereas  DLM  would  not 
do  so. 


Fig  3.23  Effects  of  Sweptback  Angle  on  and  (Af  =  0.8,  k  =  0.0,  Xo  -  0.5c) 

70  °  Delta  Wins 

A  thorough  numerical  study  has  been  conducted  on  a  70°  delta  wing  in  pitching  about  its 
midchord  at  A/  =  0.8.  The  imsteady  pressures,  forces  and  moments  are  computed  at  various 
spanwise  locations  based  on  DLM  and  ZONA6.  Figs  3.24,  3.25  and  3.26  present  the  pressure 
distribution  on  two  spanwise  stations  at  15%  half  span  (station  2)  and  95%  half  span  (station  10) 
at  reduced  frequencies  k  =  0.0,  0.5  and  10.0,  respectively.  Two  paneling  schemes  are  used,  one 
with  10  X  10  cut,  the  other  40  x  10  cut.  Several  observations  on  the  computed  pressures  in  these 
figures  are  in  order: 
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•  All  Z0NA6  results  are  valid  showing  consistent  improvement  in  pressure  details  as 
the  panel  number  increases. 

•  Almost  all  DLM  results  are  invalid  except  at  station  2  with  a  10  x  10  cut, 
corresponding  to  a  panel  aspect  ratio  of  0.43.  It  should  be  noted  that  the  panel  aspect 
ratio  for  all  other  cases  are  greater  than  unity. 

•  The  DLM  pressures  tend  to  converge  to  the  ZONA6  pressures  toward  the  aft  wing 
portions  at  station  2.  The  reason  for  this  is  not  clear.  Perhaps  the  skewness  of  the 
upstream  panels  might  also  attribute  to  the  solution  breakdown. 

An  obvious  remedy  to  use  DLM  is  to  adopt  a  different  paneling  scheme  from  the  present  one. 
This  then  would  amount  to  establishing  another  set  of  numerical  restrictions.  But  all  these  could 
be  avoided  if  one  uses  ZONA6  instead. 


10x10 


Station 


40x10 


*/c 


x/e 


Fig  3.24  In-Phase  Pressures  on  Two  Spanwise  Stations: 
70®  Delta  Wing  (M  =  0.8,  k  =  0.0,  Xo  =  0.5c) 
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Staton  (D 


Fig  3.25  Out-of  Phase  Pressures  on  Two  Spanwise  Stations: 
70“  Delta  Wing  (M  =  0.8,  k  =  0.5,  Xo  =  0.5c) 


Station 

10x10  40x10 


Fig  3.26  Out-of  Phase  Pressures  on  Two  Spanwise  Stations: 
70“  Delta  Wing  (M = 0.8,  k  =  10.0,  x„  =  0.5c 


Lastly,  Fig  3.27  presents  the  effect  of  the  panel  number  on  the  generalized  aerodynamic  forces 
and  .  It  is  seen  that  ZONA6  results  consistently  approach  to  a  converged  solution  as 
the  chordwise  panel  number  N^.  increases,  whereas  the  results  of  DLM  fail  to  do  so. 


0  TO  20  30  40  SO  0  10  20  30  40  50 


Number  of  Chordwise  Cuts  (Nx)  Number  of  Chordwise  Cuts  (Nx) 

Fig  3.27  Effects  of  Panel  Numbers  on  Cl^  and 
70“  Flat  Delta  Wing  {M  =  0.8,  k  =  10.0,  ;Co  =  0.5c) 

3.2  Integral  Equations  of  the  Linearized  Small  Disturbance  Equation 

ZONA6  and  ZONA7  solve  the  linearized  small  disturbance  equation  that  reads: 

(l  -  Mj)  4>„  +  <6,,  +  -  2  Mj  4>„  -  Mj  <»,  =  0  (3.1) 

where; 

is  the  ffeestream  Mach  number 

and; 

O  =  is  the  total  velocity  potential.  (3.2) 

Eq  3.2  shows  that  the  total  potential  O  consists  of  two  parts;  the  steady  potential  as  well 
as  the  unsteady  potential  where  is  of  the  order  of  wing  thickness  and  is  of  the  order  of 

oscillatory  amplitude  due  to  the  structural  modes.  For  stability  analysis  the  amplitude  of  the 
structural  oscillation  is  assumed  to  be  smaller  than  wing  thickness  or  at  all  times,  i.e. 

^  (3.3) 
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Substituting  Eq  3.2  into  Eq  3. 1  and  collecting  the  like-order  terms  yield  the  equations  for  and 

k- 


(l  -  mJ-)  («<„  +«*„„+  =  0  (3.4) 

(l  -  MJ)  «>,„  +  «*i„  +  (<1=  -  -  2  mJ  ^  0  (3.5) 

Eqs  3.4  and  3.5  are  the  steady  and  unsteady  linearized  small  disturbance  equations  respectively. 
The  steady  and  unsteady  potentials  are  decoupled  in  these  equations  because  Eq  3.1  is  a  linear 
equation.  Solving  the  steady  and  unsteady  flows  requires  Eqs  3.4  and  3.5  to  be  subjected  to  their 
respective  steady  and  unsteady  boundary  conditions.  The  solution  techniques  for  solving  the 
time-domain  unsteady  linearized  small  disturbance  equation  (Eq  3.5)  have  been  presented  by 
Morino  and  Tseng  (Ref  3.25)  for  subsonic  flow  and  Freedman  and  Tseng  (Ref  3.26)  for 
supersonic  flow  by  a  time  domain  Green’s  function  approach.  However,  such  a  time  domain 
Green’s  function  approach  is  complex  and  costly  in  computation.  Also,  insufficient  unsteady 
results  were  presented  in  their  work  for  proper  assessment  of  the  applicability  to  aeroelastic 
analysis. 

Consider  an  aircraft  of  interest  performs  a  simple  harmonic  motion,  the  unsteady  potential  ^ 
generated  can  be  expressed  in  a  reduced  form: 


=  (j)  e"“' 

(3.6) 

where: 

(O 

is  oscillation  frequency 

<!> 

is  the  reduced  frequency-domain  potential. 

Let: 

x' 

li 

y'  =  L  y  ,  z'  =  L  z 

(3.7) 

where; 

L 

is  the  reference  length 

(3.8) 

and: 

=  J  i-AfJ 

(3.9) 

Introducing  the  so-called  modified  potential  ^  such  that: 


<l>  =  <l> 


iXM,x' 

e 


(3.10) 


where; 


35 


II 

^8^ 

is  the  so-called  compressible  reduced  fi^equency 

(3.11) 

and: 

k  -  ^  ^ 

V 

is  the  reduced  frequency. 

(3.12) 

One  can  transform  Eq  3.5  into: 

-  2  - 

(3.13) 

^x’x'  +  ^yy 

+  n  +  A  =  0 

where: 

M  =  1 

for  <  1 

(3.14) 

//  =  -1 

for  A/oo  >  1 

(3.15) 

Applying  Green’s  theorem  to  Eq  3.13,  an  integral  solution  can  be  obtained  in  terms  of  the 
imsteady  source  and  doublet  singularity  distributions  a  and  over  the  surface  5  of  the 
configuration  of  interest.  Transforming  this  integral  solution  back  to  an  expression  for  if>,  one 
obtains: 


=  -  -J-  JJ y> 

t  n  ^ 

+  -^  JJ  %  z)  KdS 

tn  on 

where: 

■^o>  3^0’  points  (to  be  influenced) 

X,  y,  z  are  the  sending  points  (from  the  sources  and  doublets) 


£  =  4 

for  <  1 

£  =  2 

for  >  1 

a 

is  the  unsteady  source  singularity  distribution 

is  the  unsteady  doublet  singularity  distribution 

-^  =  n  •  V 
dn 

,  n  is  the  out-normal  vector 

K 

is  the  Kernel  fimction 

(3.16) 
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•  ikR 


tl 

for 

M  <  \  (Subsonic  Kernel) 

^  cos(XR) 

R 

for 

M  >  \  (Supersonic  Kernel) 

R  ^  44"  +  Rn" 

+  pi  c 

(3.17) 


(3.18) 


#  = 


r 


X  -  X 


0 


A 


V 


y  -  yp 

L 


(3.19) 


The  integral  associated  with  the  unsteady  source  singularity  is  defined  over  the  surface  of  the 
configuration  S,  where  S  can  be  expressed  as: 


Six,  y,  z,  r)  =  0 


(3.20) 


whereas  the  unsteady  doublet  singularity  over  the  configuration  and  the  wake  surfaces  S  +  IF. 
The  surface  of  a  ■typical  wing  section  and  its  associated  wake  surface  are  shown  in  Fig  3.28. 


n 


S(x,y,z,t)  =  0 


Wake 


00 


Fig  3.28  Surface  definition  of  configuration  and  wake 

The  vector  n  shown  in  Fig  3.28  represents  the  out-normed  vector  of  S  and  can  be  expressed 
as: 


n 


n^i  +  riy  j  +  n^k 


(3.21) 


where  n^,  %  and  rh  are  the  components  of  n  along  x,  y  and  z  directions,  respectively. 

For  solving  the  integral  equation  in  Eq  3.16,  ZONA6  and  ZONA7  adopt  the  so-called  panel 
method  whose  formulation  will  be  discussed  in  a  later  section. 
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3.3  Unsteady  Boundary  Condition  and  Pressure  Coefficients 

The  inviscid  fluid  flow  boundary  condition  requires  the  flow  to  be  tangential  only  to  a  moving 
body  at  all  times,  i.e.: 

^  +  F  •  V5  =  0  at  S(x.y.z.  /)  =  0  (3.22) 

dt 

where  S  is  the  body  surface  defined  in  a  wind-fixed  coordinate  system  and  V  is  the  flow 
velocity. 

Eq  3.22  is  the  generic  equation  for  all  tangency  conditions  desired  for  elastic  wing  planforms  of 
unsteady  motions.  For  elastic  bodies  performing  unsteady  motion,  however,  Eq  3.22  is  found  to 
be  unsuitable.  Attempts  in  the  past  (Ref  3.20)  to  derive  the  body  tangency  condition  using  Eq 
3.22  in  fact  failed  to  produce  useful  results.  By  useful  results,  it  is  meant  to  arrive  a  set  of 
second  order  terms  to  account  for  the  body  thickness  effect  in  addition  to  the  slender  body 
boundary  condition.  Benneker,  Roos  and  Zwaan  (Ref  3.20)  has  derived  such  a  wind-fixed 
boundary  condition.  But  their  boundary  condition  contains  several  second-order  singular  terms, 
all  associated  with  the  second  derivative  of  the  mean-flow  potential.  It  appears  that  no  clear 
resolution  was  stated  in  Ref  3.20  as  to  the  remedy  of  the  singular  terms. 

In  their  extensive  analysis,  Garcia-Fogeda  and  Liu  (Ref  3.27)  suggested  an  approach  to  adopt  the 
body-fixed  coordinate  system  that  could  totally  circumvent  the  singularity  problem.  The  flow 
tangency  condition  under  this  system  reads: 

(f  -  Fb  )  •  V5  =  0  at  S  =  0  (3.23a) 

where  is  the  velocity  due  to  body  surface  motion.  This  approach  results  in  all  regular 
second-order  terms.  The  usefulness  of  the  boundary  condition  and  the  resulting  Cp  expression 
(will  be  shown  later),  have  been  ascertained  and  adopted  by  much  of  ZONA’s  subsequent  works 
in  unsteady  subsonic  and  supersonic  aerodynamics  for  bodies  and  wing-body  configurations 
(Refs  3.28, 3.29). 

Garcia-Fogeda  and  Liu  started  from  Eq  3.23a  and  obtained  the  following  equation  for  an 
arbitrary  elastic  body  in  harmonic  motion  with  given  modes,  i.e.: 

•  «  =  Fb  (3.23b) 

where  =  (w,v,w)  is  the  unsteady  perturbation  velocity  vector  due  to  the  harmonic  motion, 
and: 


=  +  +ik\)+  n^(h^  +ikh^)  (3.24) 

is  the  so-called  “downwash  function”  on  arbitrary  bodies. 
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where: 

h^,  hy  and  fh  are  the  components  of  the  modes  with  respect  to  the  reference 
length  L  along  x,  y  and  z  directions,  respectively 

Uo  is  the  steady  perturbation  velocity  components  in  the  x  direction 


is  the  differential  operator  with  respect  to 


For  a  flat  plate  type  lifting  surface,  i.e.  Wx  =  0,  Eq  3.24  is  reduced  to  the  familiar  boundary 
condition: 

V<p  =  at  5  =  0  (3.25) 


where: 

¥^=hj+ikh^  (3.26) 

is  the  so-called  “downwash  function”  on  flat-plate  type  of  lifting  surface, 

and: 

hw  =  riyhy  +  (3.27) 

Comparing  Fw  to  Fb,  one  can  see  that  indeed  Eq  3.25  is  a  special  case  of  Eq  3.23.  Also,  while 
the  unsteady  boundary  condition  of  lifting  surfaces  is  totally  uncoupled  from  the  steady  flow 
influence,  the  unsteady  boundary  condition  of  arbitrary  body  includes  the  steady  (mean)  flow 
component  Wo-  Computing  the  steady  flow  components  requires  the  solution  of  the  steady  small 
disturbance  equation  (Eq  3.4)  with  the  well-known  steady  boundary  condition: 

•  n  = -(F  •  «)  (3.28) 


where  V  is  the  free-stream  velocity  vector. 

Based  on  the  work  of  Garcia-Fogeda  and  Liu  (Ref  3.27),  the  unsteady  pressure  coefScient  on  the 
body  expressed  in  the  body-fixed  coordinate  reads: 


Cp  =  -2  C„ 


+  u^+hy  «o)  +  +  w,,  w 

+  +  hy  v„  +  ik{<f>  +  /ix  Wo  +  hy  v„  +  K 


where: 


Co 


1 


(lu^+ul+vl^-wl 


(3.29) 


(3.30) 
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Uo,  Vo  and  Wo  are  the  steady  perturbation  velocity  components 

^,u,v  and  w  are  the  unsteady  perturbation  potential  and  velocity  components 

y  is  the  specific  heat  ratio 

For  a  flat-plate  type  lifting  surface  at  zero  angle  of  attack,  the  steady  perturbation  velocity 
components  are  zeros,  i.e.  Wo  =  Vo  =  Wo  =  0.  Eq  3.29  is  the  reduced  to  the  well-known 
pressure  coefficient  expression  of  lifting  surfaces: 

Cp  =  -2  (w  +  ik^)  (3.31) 

However,  the  conventional  pressure  coefficient  of  lifting  surface  is  an  expression  in  terms  of  the 
difference  of  pressure  between  its  lower  and  upper  surfaces.  Expressing  Eq  3.31  in  terms  of 
pressure  difference  on  the  mean  surface  gives: 

ACp  =  -2  (Am  +  ikA^)  =  "2  (3.32) 

where  A(-)  denotes  the  difference  between  the  lower  and  upper  surfaces  of  (•) . 

It  should  be  remarked  that  the  xmsteady  pressure  Cp  for  arbitrary  bodies  Eq  3.29  involves 
coupling  terms  with  the  perturbation  velocities  of  the  steady  mean  flow,  which  brings  in  the 
thickness  effects.  By  contrast,  the  unsteady  pressure  ACp  for  lifting  surface  (Eq  3.32)  is 

imcoupled  fi'om  the  steady  mean  flow  terms. 

3.4  Paneling  Scheme  for  Aircraft  Configurations 

The  surface  of  an  aircraft  configuration  is  broadly  divided  into  two  categories:  the  body-like 
components  and  wing-like  components  (see  Fig  3.29).  The  body-like  components  (BLC) 
consists  of  fuselage,  external  stores,  tiptanks,  etc.  Each  BLC  is  divided  into  NG  number  of 
segments  by  x  =  constant  cuts.  Each  segments  consists  of  NR  x  NQ  number  of  surface  boxes,  all 
of  which  are  flat  panel  elements.  That  is,  each  segment  is  subdivided  into  NR  rings  by  NR+\ 
cross  sections,  obtained  by  cutting  the  segment  with  x  =  constant  planes.  Next,  a  circumferential 
(radial)  cut  is  applied  to  this  same  segment  yielding  NQ  boxes  on  each  ring.  There  are  NQ+\ 
circumferential  points  to  define  the  NQ  boxes.  The  order  of  these  NQ^\  points  begins  from  the 
bottom  meridian  coxmting  in  the  clockwise  direction  looking  downstream.  Let  NB  represent  the 
number  of  boxes  of  the  total  body-like  components  for  an  aircraft,  and  let  NBODY  denote  the 
number  of  body-like  components.  Thus,  NB  can  be  expressed  as: 

NBODY  NG, 

ffB  =  Y,  (3.33) 

/=1  /=1 
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Chordwise 


Fig  3.29  Aircraft  Components  Showing  Body  and  Wing  Segment/Box 

The  wing-like  components  (WLC)  consists  of  thin  surfaces  whose  spanwise  cross  section  can  be 
represented  by  airfoil-like  thickness  distribution  such  as  wings,  canard,  fins,  pylons,  la.unchers, 
etc.  A  wing  component  is  defined  by  a  lifting  surface  enclosed  by  a  straight  leading,  trailing  and 
two  side  edges.  Each  wing  component  is  then  divided  into  NS  streamwise  strips  by  NS+l  liiies 
parallel  to  the  x-axis  called  “spanwise  divisions”,  and  each  strip  is  subdivided  into  NC  chordwise 
trapezoidal  boxes  by  NC+l  spanwise  lines  called  “chordwise  divisions”.  Let  iVTf^be  the  number 
of  boxes  of  the  total  WLC  of  an  aircraft,  and  let  NWING  denote  the  number  of  wing  segments. 
Then  NW  can  be  expressed  as; 

miffG 

NW  =  Y^NS.xNC,  (3-34) 

It  should  be  noted  that  for  unsteady  linear  aerodynamics  the  thickness  effects  of  the  thin  lifting 
surfaces  are  of  first  order.  This  is  to  say  that  for  the  wing-like  components  there  is  no 
requirement  of  modeling  the  thickness  distribution  on  the  lifting  surfaces.  Thus,  all  wing-like 
components  are  assumed  to  be  flat  plate  in  the  present  paneling  scheme.  Since  source  singularity 
is  usually  used  to  simulate  the  thickness  effects  whereas  doublet  singularity  to  generate  lift,  the 
flat-plate  type  of  lifting  surfaces  require  no  source  singularity  and  only  the  doublet  singularity  is 
distributed  on  the  mean-plate  of  the  wing-like  components  and  their  associated  wake  siorfaces. 
On  the  other  hand,  the  present  panel  scheme  for  bodies  takes  up  a  simplified  formulation  in 
which  only  the  source  singularity  is  distributed  on  the  body-like  components.  This  will  greatly 
simplify  the  effort  for  deriving  the  solutions  of  the  source  and  doublet  integrals.  With  this 
paneling  scheme,  Eq  3.16  becomes: 
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(3.35) 


where: 


^  JJ  l^<f>{x,y,z)  -^KdS 


(3.36) 


wing + wake 


and: 


is  the  potential  influence  due  to  the  wing-like  components 


JJ  y,  z)  KdS 

^  ^  body 


(3.37) 


is  the  potential  due  to  the  body-like  components 


The  derivations  of  and  solution  are  discussed  on  the  next  section. 

3.5  Discretization  of  the  Source  and  Doublet  Integrals 

Doublet  singularity  distribution  on  the  wing-like  components  and  their  associated  wake  surfaces 
implies  that  the  panel  modeling  is  required  not  only  for  the  wing-like  components  but  also  for  the 
waJce  surfaces.  This  will  greatly  increase  the  modeling  effort  due  to  the  large  domain  of  the 
wake  surface  that  starts  from  the  wing  trailing  edge  and  extends  to  downstream  infinity.  To 
circumvent  this  problem,  one  introduces  the  acceleration  potential  v  which  is  directly 
proportional  to  a  perturbation  pressure  field  related  to  the  potential  ^  by: 


y/  —  - —  + 

dx 


(3.38) 


Comparing  Eq  3.38  with  Eq  3.32  results: 
AC.  =  -2y/  =  -2  f 


(3.39) 


Eq  3.39  establishes  a  transformation  from  the  unsteady  doublet  singularity  A^^  to  the  pressure 
coefficient  difference  ACp ,  Applying  this  transformation  to  Eq  3.36  gives: 


(  -^O  ’  .Vo  »  ^0  ) 


JJ  K(4,^,OdS 


(3.40) 


where: 


(3.41) 
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and: 

X  =:  -  oo  for  A/  <  1 

X  =  fox  M>  1 

Since  the  wake  surface  cannot  sustain  force,  i.e.  ACp  =  0  on  wake  surface,  the  integral  equation 
of  (f>^  (Eq  3.40)  is  applied  on  the  mean-plane  of  the  wing  components  only.  The  integral  on 
wake  surface  vanishes  automatically.  Therefore,  with  the  acceleration  potential  transformation, 
the  integral  equation  of  expressed  directly  in  terms  of  ACp  and  requires  no  panel 

modeling  of  the  wake  surface. 

Based  on  the  paneling  scheme  depicted  in  Fig  3.29,  the  integrals  of  Eqs  3.37  and  3.40  can  be 
discretized  into  many  elementary  Kernel  integrals  associated  with  each  aerodjmamic  box.  On 
each  aerodynamic  box,  a  constant  singularity  distributed  is  assumed.  For  the  /  control  point 
located  either  on  the  wing  or  body  box,  its  perturbation  potential  ^  can  be  approximated  by: 


+  Z  AC, 


j^\  j^l 

(  i  =  1,2,. ..,NB+NW) 


pj 


(3.42) 


where: 


=  .  _i_  ff  K  dn 

■y  Eji 


(3.43) 


is  the  elementary  source  Kernel  integral  representing  the  potential  influence 
coefficient  at  the  control  point  due  to  they*  body  box 


Jj  e-'*^^  K  drj  d^ 


(3.44) 


is  the  elementary  pressure  Kernel  integral  representing  the  potential  influence 
coefficient  at  the  i*  control  point  due  the y*  wing  box 


<7j  is  the  unknown  soiuce  strength  of  the  y  body  box 
ACp^  is  the  unknown  pressure  jump  on  they*  wing  box 

and  Atjy  represent  the  boundary  of  the  y*  box  relative  to  the  /*  control  point 
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Likewise,  the  unsteady  perturbation  velocities  w/,  v„  w,,  can  be  expressed  as; 


NB  Ntr 

^  (3-45) 

j^\ 

(  i  =  1,2,...,NB+NIV) 

where: 

(“b,.vb.,«'b.  )  =  —  n  ^ 

and: 

(t'w..v„,,«'w.)  =  7^-  JJ  v{e'‘'’*^}dvdf  (3.47) 

are  the  velocity  influence  coefficients  at  the  control  point  due  to  they*^ 
wing  box 

The  potential  and  velocity  influence  coefficients  of  the  body  box  and  the  wing  box  are  obtained 
by  solving  the  elementary  kernel  integrals  (Eq  3.44)  and  the  elementary  pressure  kernel  integrals 
(Eq  3.45)  respectively.  The  detailed  derivations  of  the  solutions  of  these  integrals  can  be  found 
inRefs  3.2, 3.15  and  3.19. 

3.6  Matrix  Equations  for  the  Solution  of  Unsteady  Pressure 

Applying  Eqs  3.42  and  3.43  at  number  of  control  points,  one  can  construct  the  influence 

coefficient  matrices  P/C,  UIC,  VIC  and  WIC  as: 

m  =  [«c] 

=  '“"^5  {acJ 
w  =  ivin 

(w)  =  mc\  T"  I  (3.48) 

where  ,  {u} ,  {v}  and  {w}  are  the  perturbation  potential  and  velocities  of  all  control  points, 
and: 
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is  the  potential  influence  coefficient  matrix. 


[PIC]  = 

(?^b)bb 

(^w)wB 

(^b)bw 

(^w)ww  . 

[UIC]  = 

(^b)bb 

(*^w)wB 

.  (%)bw 

(“w)ww  _ 

[WIC]  = 

(^b)bb 

(^w)wB 

_  (^b)bw 

(^w)ww 

[VIC]  = 


('’b)bb 

(^b)bw 


(''w)wB 

(^w)ww 


are  the  velocity  influence  coefficient  matrices. 


(3.49) 


().B 

(  )bw 

(  )wB 
(  )ww 


is  the  influence  at  the  body  control  points  due  to  the  body  boxes 
is  the  influence  at  the  wing  control  points  due  to  the  body  boxes 
is  the  influence  at  the  body  control  points  due  to  the  wing  boxes 
is  the  influence  at  the  wing  control  points  due  to  the  wing  boxes 


Applying  the  boundary  conditions  of  the  wing  boxes  (Eq  3.27)  and  the  body  boxes  (Eq  3.24) 
gives; 


[NIC] 


(3.49a) 


where: 

[MC]  =  i?  •  « 

=  [nJ[WC]  +  [n^WIC]  +  [n,][WIC] 

is  an  {NB+NW)  by  {NB+NW)  square  matrix  and  is  named 
Normal  Velocity  Influence  Coefficient  Matrix 

u  =  (m,v,w) 

is  unsteady  velocities  on  wing-body 

Inverting  [NIC]  gives  the  singularity  strengths  of  the  body  and  the  wing  boxes,  i.e.: 


For  a  given  set  of  downwash  functions  Fw  and  Fb,  unknown  strengths  o"s  and  ACp ’s  can  be 
solved  from  Eq  3.50.  Once  ct’s  and  ACp ’s  are  known,  the  perturbation  potential  and  velocities 
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,  {u} ,  {v}  and  {w}  can  be  evaluated  according  to  Eq  3.48.  Consequently,  the  unsteady 
pressure  coefficient  on  body  boxes  can  be  calculated  according  to  Eq  3.29. 

Notice  that  the  boundary  condition  and  pressure  coefficient  of  the  body  boxes  involve  the  steady 
perturbation  velocities  «o,  Vo  and  Wo.  To  evaluate  them  form  Eq  3.50,  it  is  required  to  solve  the 
steady  aerodynamics  of  the  body.  Now,  the  unsteady  Kernel  functions  of  Eq  3.17  can  be  readily 
reduced  to  its  counterpart  in  the  zero  frequency  limit.  Thus,  the  influence  coefficient  matrices  of 
the  steady  aerodynamics  can  adopt  directly  from  these  of  the  unsteady  aerodynamics  by  setting 
k  =  0.  The  steady  perturbation  velocities  can  then  be  obtained  by  applying  the  steady  boundary 
condition  presented  in  Eq  3.28. 

Fig  3.30  presents  the  solution  procedure  in  obtaining  the  unsteady  pressure  coefficients  on  wing 
and  body  boxes.  This  procedure  is  proceeded  as  follows; 

-  Solve  for  steady  perturbation  velocities  Uo,  v©  and  Wo. 

-  Construct  Influence  Coefficient  matrices  {UIC\,  [VIC\,  [WIC\,  {PIC\  and  [NIQ  according 
to  Eq  3.49.  Note  that  this  step  takes  up  most  of  Ae  computation  time  in 
ZONA6/ZONA7.  However,  being  independent  of  the  mode  shape,  they  can  be  obtained 
and  saved  once  and  for  all  for  a  given  A/ and  k  pair. 

-  Construct  Downwash  functions  and  of  Eqs  3.24  and  3.26;  perform  matrix 
decomposition  on  [MC]  and  solve  for  crand  ACp  from  Eq  3.50. 


-  Compute  unsteady  velocities  w,  v  and  w  and  unsteady  potential  ^  according  to  Eq  3.48. 

-  Compute  imsteady  pressure  Cp  on  body  according  to  Eq  3.29. 
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Fig  3.30  Flow  Chart  of  Computation  and  Calculation  Procedure 
for  Unsteady  Pressures 


3.7  Construction  of  Aerodynamic  Influence  Coefficient  (AIC)  Matrix 

The  flow  chart  shown  in  Fig  3.30  gives  a  step-by-step  computational  procedure  for  unsteady 
pressure  calculation,  but  it  does  not  provide  the  AIC  matrix  to  directly  relate  the  given  structural 
mode  shape  to  the  unsteady  aerodynamic  forces.  To  construct  an  AIC  matrix  for  the  wing-like 
components  is  rather  straightforward.  The  force  acting  normal  to  the  wing  boxes  can  be  obtained 
by  multiplying  the  area  of  the  wing  boxes  to  the  unsteady  pressure  on  wing  boxes.  This  gives: 

(W  =  ?.[«  A,][MCr'  {J»| 


(3.51) 


where; 


^  ^  L[^^C3bw  [njc]^ 

is  the  normal  velocity  influence  coefficient  matrix 
g„  is  the  dynamic  pressure 

Ay,  is  a  diagonal  matrix  containing  the  area  of  the  wing  boxes 
is  the  normal  force  vector  on  wing  boxes 


The  normal  forces  on  the  body-like  components  are  more  complicated  that  that  of  the  wing-like 
components.  Rewriting  Eq  3.29  in  matrix  form  and  multiplying  the  body  box  area  give; 


[Ab)  {C,} 

■=  9.  [[Bbb]  DBwb]]  Cj  +  {d) 

I'wJ 


(3.52) 


where; 

A  B  is  a  diagonal  matrix  containing  the  area  of  the  body  boxes 

L  body  is  the  normal  force  vector  on  the  body  boxes 

[[Bbb]  [Bwb]]  =  [AbI  [-2CJ  [  [[WC,b)  [WC»b1] 

+  K]  [F/CbbI  [WC^l] 

+  [F/Cb,]  [WCwbI] 

+  [<■*]  [[WCbbI  [WC^lj] 

where  d  contains  the  terms  in  Eq  3.29  associated  with  the  given  mode  shapes, 
{d}  =  -2  Ab  Co  [  (1  +  Mo )  (K  "o  ^o)  +  K  \ 

+  ik(Ji^  Mo  +  Ay  Vo  +  Wo)  ] 


(3.53) 


(3.54) 


Combining  Eqs  3.51  and  3.54  gives; 


q.  [NIC] 


in 

iFwJ  loj 


(3.55) 


where; 
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fel  =  r  [®“|  1“™]]  (3.56) 

L[*]  KlJ 

Eq  3.55  establishes  a  relationship  between  the  boundary  condition  of  box  and  wing  components 
(Eqs  3.24  and  3.25)  to  the  normal  forces  on  the  box  and  wing  boxes.  However,  it  is  desirable  to 
express  Eq  3.55  in  the  standard  format  of  AJC  matrix  equation  as; 

{LO  =  9.  {AIC]  {h}  (3.57) 

where  h  and  Lh  represents  respectively  expression  for  the  mode  shapes  and  force  vectors  on  the 
body  and  wing  boxes. 

To  derive  the  Eq  3.57  from  Eq  3.55,  it  is  required  to  introduce  the  J-set  and  K-set  degrees  of 
freedom  of  the  aerodynamic  boxes.  This  will  be  discussed  next. 

3.8  J-Set  and  K-Set  Aerodynamic  Degrees  of  Freedom  for  AIC  Matrix 

The  J-set  aerodynamic  degrees  of  freedom  (d.o.f.)  is  simply  the  number  of  aerodynamic  boxes. 
For  instance,  the  size  of  vectors  and  matrices  inEq  3.55  are  all  in  the  oxdsxoi  J-set  d.o.f ,  i.e.: 

J-set  d.o.f.  =  NB  +  NW  (3-58) 

By  examining  the  boundary  conditions  of  the  body  and  wing  boxes,  it  can  be  realized  that  mode 
shapes  h  in  Eq  3.57  contains  six  components  on  each  aerodynamic  box,  namely 
h^,hy,h^,hl^,  h'y  and  h[ .  This  implies  that  the  size  of  h  is  in  the  order  of  6  x  J-set.  Defining  K- 

set  =  6  X  J-set  as; 

K-set  =  6  X  J-set  =  6  x  (NB+NW)  (3-59) 


gives  the  mode  shapes  h  as; 


(3.60) 


where  i  denotes  the d.o.f, 
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and  (•)'  represents  the  derivative  of  (•)  with  respect  to 


Rewriting  the  vectors  {Fg}  and  {d}  in  Eq  3.55  in  matrix  forms,  one  can  redefine  Eq  3.55  as: 


q.  [HiMC]-' 


(3.61) 


where  F  and  D  are  the  K-set  by  J-set  complex  matrices  containing  the  steady  velocity 
components  «o,  Vo  and  Wo  and  the  normal  vector  components  «x,  %  and  riz  as  well  as  the  reduced 
frequency  k. 


In  general,  the  force  at  each  grid  point  of  the  structural  finite  element  method  has  six 
components;  forces  and  moments  along  the  x,  y  and  z  directions.  In  order  to  be  compatible  with 
the  structural  finite  element  method,  it  is  required  to  convert  the  normal  force  vector  in  Eq  3.61 
to  the  three  force  components  (F^,Fy,F^)  and  the  three  moment  components 
along  the  x,  y  and  z  directions.  This  can  be  done  by: 


(3.62) 


Combining  Eqs  3.61  and  3.62  gives: 

{L,}  =  [AlC]  {h}  (3.63) 

Finally,  we  have  derived  a  K-set  by  K-set  square  AlC  matrix  that  directly  relates  the  structural 
mode  shapes  to  the  aerodynamic  forces.  It  should  be  noted  that  the  structural  mode  shape  (h)  in 
Eq  3.62  is  defined  at  the  aerodynamic  boxes.  Transferring  the  structural  mode  shapes  from  the 
structural  finite  element  grid  points  to  the  aerodynamic  boxes  requires  the  generation  of  a  spline 
matrix.  This  will  be  discussed  in  Section  6. 


3.9  Body  Wake  Effect 

The  wake  generated  by  flow  separation  at  the  tail  section  or  behind  the  base  of  a  body-like 
component  has  considerable  influence  to  the  wing-body  or  body-alone  aerodynamics.  ZAERO 
provides  a  non-iterative  approach  for  the  body  wake  modeling.  Fig  3.31  presents  the 
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comparisons  of  the  ZONA6  results  with  and  without  the  wake  modeling  with  the  experimental 
pressure  distributions  along  the  blimt-base  body  in  incompressible  flow. 


Fig  3.31  Comparison  of  Surface  Pressure  Distribution  for  a  Blunt  Body  (Ls/d  -  5)  at  a  -  0®,  M  - 

0  and  Base  Pressure  =  -  0.169 

It  is  seen  that  the  present  (no  wake)  result  begins  to  deviate  from  the  measured  data  at  a  forebody 
portion  of  0.3  Ib,  where  Lb  is  the  body  length.  By  contrast,  the  present  (with  wake)  result  shows 
remarkably  close  agreement  with  measured  data.  This  comparison  clearly  demonstrates  the 
importance  of  the  body  wake  modeling  for  accurate  aerodynamic  prediction. 

This  non-iterative  body  wake  modeling  requires  a  given  base  pressure  coefficient  a  priori. 
Therefore,  the  pressure  of  the  adjacent  boxes  to  the  body  base  is  given  by  Cp  =  .  Imbedded 

singularities  are  placed  in  the  assigned  proximity  of  the  body  base  regime  to  simulate  the  exterior 
wake  flow.  For  steady  flow  the  imbedded  singularity  is  a  source,  whereas  for  unsteady  flow  it  is 
a  doublet.  Next,  the  constant  pressure  condition  is  imposed  at  the  body  base.  This  boundary 
condition  can  be  expressed  as: 

C  =C  and  Cp=0  at  S(x=L^,^z)  =  0  (3.64) 

Po  Pb#«  r 

where  C„  and  C.  are  the  steady  and  unsteady  pressure,  respectively. 

Vq  * 
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The  combination  of  this  constant  pressure  condition  and  the  steady  and  unsteady  boundary 
conditions  yields  a  compact  expression  in  terms  of  the  velocities,  the  mode  shape  and  its 
derivative,  i.e.: 

a  (f)  -  b  (3.65) 

where  a  and  h  are  both  scalar  operators. 

a  =  Vo  '  ^ 

^  .  V  +  •  V 

dx 

c  is  a  function  of  the  mode  shape 

and  ^  are  the  steady  and  unsteady  potential  due  to  the  imbedded  singularity, 
respectively. 

The  location  of  the  imbedded  source  and  the  imbedded  doublet  are  placed  along  the  extended 
body  axis  at  x  =  and  x  =  +  Qi,  respectively,  in  the  region  near  the  base.  The 

values  of  and  qi  are,  in  general,  related  to  the  wake  length  which  is  yet  to  be  solved  as  part 
of  the  solution,  ordinarily,  to  be  solved  through  an  iteration  procedure.  However,  our  experience 
showed  that  empirical  guidelines  can  be  set  up  in  which  q^  can  be  confined  to  a  width  of  0.2  ~ 
0.4  Lb  and  q^  ofo.05~0.15  Lb  for  any  circular  blunt  base  (Z-b  being  the  body  length).  By 
numerical  experiment,  it  is  found  that  the  converged  solutions  are  rather  insensitive  to  the  precise 
location  q^  and  q^  as  long  as  they  are  placed  within  the  width  given.  Fig  3.32  shows  the  side 
views  of  the  computed  base-flow  of  a  blunt  body,  according  to  the  present  no-wake  and  the 
present  with-wake  model.  The  direction  and  the  length  of  the  arrows  indicate  the  flow  direction 
and  the  velocity  magnitude,  respectively.  As  expected,  the  no-wake  flowfield  velocities  appear 
to  be  nearly  uniform  showing  their  immediate  return  to  the  freestream  condition.  It  can  be  seen 
that  the  present  wake  model  appears  to  closely  simulate  the  wake  profile.  The  dashed  line  is 
indeed  the  dividing  streamline  which  defines  the  wake  closme,  or  the  computed  wake  shape.  It 
is  clearly  shown  that  the  stagnation  point  at  the  end  of  the  wake  is  automatically  captured  by  the 
present  calculation.  The  computed  wake  flow  is  only  physically  meaningful  in  the  outer  wake 
flow  region,  whereas,  the  computed  flow  inside  the  wake  is  of  no  physical  significance. 
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a)  without  wake  model 


Fig  3.32  Computer  Wake  Shape  for  a  Blunt  Body  at  M  =  0  and  a  =  0“; 

a  Meridian-Plane  View 

With  the  boundary  condition  expressed  in  Eq  3.65  for  body  wake  effects,  the  unknown  strength 
of  the  source  singularity  ct  on  body  boxes,  ACp  on  wing  boxes,  and  the  imbedded  wake 

singularity  can  be  solved  by  the  following  MC  matrix  equation: 

■  [^CJbb  [MCIbw  [MC]b  -wake  ^  1  [  1 

[MCI™  [NIC]^  [NIC]^.^  I  AC,  =  I  F„  I  (3.66) 

„  (^^^)wakc-B  wake-W  (^^)»«ke-wak«  J  t /^wake  J  L  ^  J 

where: 

[MC]B.^e  and  [^C]w.,«.ke  represent  the  NIC  induced  on  the  body  and  wing  boxes  by 
the  wake  singularity. 

(«^^)wake.w  and  are  the  constant  pressuTe  condition  imposed  on 

the  body  boxes,  wing  boxes  and  the  wake  singularity,  respectively. 

Once  the  unknowns  a ,  ACp  and  are  solved  from  Eq  3.66,  the  unsteady  potential  ^  and 
velocities  u,  v  and  w  can  be  obtained  from  Eq  3.48. 

Fig  3.33  presents  the  unsteady  pressures  along  the  tip-tank  of  the  NLR  wing-tip-tank-pylon-store 
configuration  (Fig  3.34)  at  an  azimuthal  angle  6= 202.5° .  The  body  wake  effects  on  the  wing- 
body  configuration  can  be  clearly  seen  in  the  ZONA6  results  using  the  with-  and  without-wake 
models. 


Fig  3.33  Unsteady  Pressure  Distribution  Along  the  Tip-Tank  of  NLR  Wing-Tip-Tank 
Configuration;  M  —  0.45,  a  =  0°,  A  =  0.305,  =  0.15  Cr,  and  0  =  202.5® 


AERODYNAMIC  MODEUNG 
WING  BOXES  *90 
PYLON  BOXES -24 
TANK  PANELS -264 
STORE  PANELS -216 


Fig  3.34  NLR  Wing-Tip-Tank-Pylon-Store  Configuration  Showing  Paneling  Scheme 


54 


It  is  seen  that  the  present  no-wake  in- 
phase  Cp  appears  to  deteriorate  towards 
the  tail  of  the  tip-tank,  whereas,  that  with 
wake  appears  to  correlate  best  with  the 
NLR  measurements.  This  is  expected, 
because  the  flow  is  known  to  separate  at 
the  rear  in  this  case,  according  to  Ref 
3.20.  For  the  out-of-phase  Cp,  the  NLR 
computed  results  (Ref  3.20)  show  a 
discrepancy  with  the  measurement 
starting  from  the  mid  body. 

The  computed  steady  pressure 
distribution  along  the  store  at  three 
azimuthal  angles,  namely, 

0  =  90,  180  and  270®,  are  presented  in 
Fig  3.35.  The  solution  with  wake  is  seen 
to  have  better  agreement  with  the 
measured  data  that  without  a  wake.  This 
leads  one  to  believe  that  the  flow  might 
have  separated  at  the  body  tail  in  this 
experiment,  although  no  such 
information  was  supplied  in  Ref  3.20. 

Fig  3.36  presents  the  imsteady  pressure 
distributions  along  the  store  at  two 
azimuthal  angles:  1)  0  =  157.5®  and  2) 
0  =  292.5° .  From  the  present  computed 
results,  it  is  seen  that  adding  the  wing-tip- 
tank-pylon  has  altered  the  imsteady 
pressure  distributions  on  the  store 
substantially.  In  all  cases,  the  present 
results  with  or  without  wake  are  in  good 
agreement  with  the  measured  data.  In 
particular,  the  solution  with  wake  seems 
to  show  better  agreement  that  the  no¬ 
wake  solution  for  the  in-phase  Cp,  and 
there  the  solution  also  shows  a  trailing- 
edge  peak  due  to  pylon-store  interaction 
(at  0  =  292.5®). 


Fig  3.35  Steady  Pressure  Distribution  Along  the  Store  of 
NLR  Wing-Tip-Tank-Pylon-Store  Configuration; 

M  =  0.45  and  a  =  0“: 
a)  0  =  90%  b)  0  =  180°  and  c)  0  =  270° 


On  the  other  hand,  NLR’s  computed  results  in  the  out-of-phase  Cp  indicates  a  solution  departure 
from  the  measured  data  starting  from  the  mid  body  section,  and  further  deteriorates  towards  the 
body  tail.  By  contrast,  it  is  seen  that  both  present  solutions  results  in  much  better  agreement  with 
the  measured  data,  especially  near  the  body  tail. 
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Fig  336  Unsteady  Pressure  Distribution  Along  the  Tip-Tank  of  NLR  Wing-Tip-Tank 
Configuration  at  a)  6  =  157.5°  and  b)  0  =  292.5°; 

M  -  0.45,  a  =  0°,  k  =  0.305  and  =  0.15  cr, 


3.10  Technique  of  Minimizing  Spurious  Waves  for  Supersonic  Body  Boxes 

High  order  panel  methods  usually  adopt  the  quadratic  singularity  distribution  on  the  body  boxes 
(e.g.  PANAIR,  Ref  3.30).  By  comparison,  the  constant  source  singularity  distribution  on  the 
ZONA6/7  body  boxes  should  be  considered  as  a  low  order  panel  approach  for  body  unsteady 
aerodynamics.  For  steady  supersonic  flow  over  a  body,  it  is  found  that  these  constant  source 
singularities  on  body  boxes  could  generate  the  nonphysical  internal  Mach  waves  or  “spurious” 
Mach  waves  inside  the  body-like  components.  These  spurious  Mach  waves  tend  to  accumulate 
as  they  propagate  dovmstream  within  a  closed  body.  As  a  result,  this  accumulation  of  waves 
often  destabilize  the  external  flow  solutions  and  hence  generate  large  error  in  the  solution  over 
the  aft  part  of  a  closed  body.  This  problem  can  be  clearly  seen  in  Fig  3.37  for  a  supersonic  flow 
over  a  cone-cylinder-cone  body.  In  contrast  to  the  exact  solution  of  Karman-Moore  (Ref  3.31), 
the  uncorrected  ZONA7  solution  (without  wave  minimization)  deteriorates  to  the  extent  that  the 
solution  at  the  aft-cone  becomes  totally  unacceptable. 

To  examine  how  the  spurious  waves  are  generated  in  supersonic  flow,  first  consider  a  point  P 
located  on  the  x-z  plane  shown  in  Fig  3.38. 
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Fig  3.37  Cone-Cylinder-Cone  at  M  =  2.0  and  a  =  0°  Using  ZONA7  without  Wave  Minimization 


z 


Fig  3.38  Opposite  Points  from  the  Panel  Having  the  Same  xjower 
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The  steady  normal  perturbation  velocity  at  P,  — from  the  body  box  (assuming  a  unit  strength 
source  singularity  is  distributed  on  the  box)  at  =  42  reads: 


5  r^«pp«  f^«pp«  1 

JjCu-,-.  jj 


&Q  J^iowcr 

l^^upper 

d  f^^UppCT  .  TJ 

I  -Sin 


dy 


dz^ 


dx 


where: 


|  =  Xo-x,  ri=^y^-y  and  C  = 


(3.67) 


and  Xupper,  J^iower,  Supper  and  >iower  denote  the  boundary  of  the  body  box. 

Eq  3.67  is  the  steady  supersonic  source  integral  of  the  body  box  obtained  by  setting  k=Q  and  p 
=  1  in  Eq  3.37.  If  the  body  box  has  a  large  width  so  that  puppet  and  yiower  are  located  forward 
outside  of  the  forward  Mach  cone  emitted  from  P,  the  integration  limits  of  Eq  3.67  is  the 
intersection  of  the  forward  Mach  cone  R  =  0  and  the  panel  z  =  0,  i.e.: 

y  -  (x^-x)  -  zl  at  a  particulars  (3.68) 


Therefore,  the  integrand  of  Eq  3.67  has  the  value  of: 
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(3.69) 


The  upper  limit  of  the  x  integral  is  at  the  panel’s  leading  edge,  i.e.  x^^^^  =  0 .  However,  based  on 
Eq  3.68  for  y  =  0,  one  can  conclude  that  -  |z„  | .  This  proves  that  for  two  points  located 

opposite  of  the  panel  with  the  same  x  location,  they  have  the  same  value  of  . 


Therefore,  Eq  3.68  is  resulted  as: 


^  I^upper-^werl 

0Z„  5z„ 


(3.70) 


=  -sign(zj7r 


Eq  3.70  shows  that  for  a  point  P  located  above  the  box  (point  Pi  in  Fig  3.38),  it  receives  a 

1 

constant  perturbation  velocity  =  -  —  whereas  a  point  P  located  below  the  box  (point  P^  in 

dz„  2 
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Fig  3.38),  it  receives 


1 

=  — .  It  also  proves  that  constant  source  panel  generates  a  jump  of 

&o  2 


across  the  box  from  Zo  =  O'  to  Zo  =  0‘*'. 


To  conclude  the  above  discussion,  consider  a  point  on  the  x  -  z  plane  travels  along  the  x 
direction.  This  point  receives  no  perturbation  from  the  box  until  its  forward  Mach  cone  hits  the 
box  leading  edge.  For  example,  at  the  canopy  sill  line  (Fig  3.39),  a  finite  jump  of  perturbation 


velocity  occurs  (i.e. 


±  sign  selection  depends  on  the  upper  or  lower  Mach  cone).  In 


other  words,  the  leading  edge  of  the  body  box  generates  two  waves;  the  outer  wave  goes  over  the 
upper  side  and  the  inner  wave  goes  under/to  the  lower  side  of  the  box.  The  outer  (exterior)  wave 
generated  from  the  leading  edge  of  the  box  propagates  into  the  exterior  flowfield  and  actually 
forms  part  of  the  required  “physical”  linearized  supersonic  exterior  flow  solution.  However,  the 
inner  (internal)  wave  generates  a  set  of  inner  wave  train  within  the  interior  domain  of  the  body. 
It  is  this  iimer  wave  train  that  creates  the  spurious  velocity  discontinuities  which  is  responsible 
for  the  solution  deterioration  at  the  aft-cone. 


ZONA6  and  ZONA7  approximate  a  body  surface  by  a  set  of  fairly  small  boxes  that  are 
approximately  contiguous  and  “flat”  in  curvature.  Clearly,  if  such  a  geometric  discretization  is 
coupled  with  discontinuous  singularity  distribution  such  as  constant  source  strength,  then  the 
amoimt  of  cancellation  of  the  spurious  waves  between  boxes  must  depend  on  the  differences 
between  their  source  strengths.  For  flat  boxes  with  constant  source,  to  keep  such  strength 
difference  small  requires  the  change  in  streamwise  geometric  slope  between  adjacent  boxes  to  be 
kept  small.  This  requirement  sets  the  size  criteria  for  box  discretization  which  should  minimize 
the  undesirable  spurious  velocity  discontinuities  to  be  generated  in  the  external  flowfields. 
Without  wave  minimization,  ZONA7  usually  works  well  only  for  bodies  with  a  simplified  body 
model,  in  which  the  actual  geometry  of  the  body  should  be  smoothed  out.  For  engine  inlets,  they 
should  be  sealed  and  smoothed  out  as  well. 

Next,  consider  a  body  with  a  realistic  geometry  comer  representing  a  canopy  sill  line  as  shown  in 
Fig  3.39.  The  canopy  sill  line  generates  a  set  of  strong  Mach  waves  hence  results  in  a  larger 
jump  in  source  strengths  between  the  two  adjacent,  dividing  the  boxes. 


Fig  3.39  Propagation  of  Spurious  Wave  generated  by  the  Real  Geometry  Comer 
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When  the  internal  wave  hits  point  P\,  the  panel  at  Pi  must  produce  a  jump  in  source  strength  in 
order  to  overcome  the  velocity  discontinuity  generated  by  the  internal  wave.  In  turn,  it  produces 
another  internal  wave  (also  another  exterior  wave)  which  travels  to  downstream  and  it  hits  point 
Pi.  As  a  result,  the  box  at  Pi  must  produce  the  same  type  of  internal  wave  that  will  hit  P3  and  so 
on.  This  “wave  train”  process  means  that  the  original  internal  wave  from  the  canopy  sill  line  will 
propagate  downstream  inside  the  body  and  the  subsequent  waves  generated  will  be  accumulated 
for  a  closed-end  body.  For  a  realistic  fuselage-like  body  with  many  geometric  comers,  numerous 
internal  waves  generated  by  them  will  be  accumulated  downstream  and  will  likely  be  further 
accentuated  if  two  of  these  waves  land  on  one  single  box.  Clearly,  these  numerous  internal 
wave-trains  often  degrade  the  (external)  flowfield  quality  to  an  extent  that  the  resulting  flowfield 
solution  becomes  unacceptable. 

On  the  basis  of  the  preceding  discussions,  it  should  now  be  clear  that  many  velocity 
discontinuities  can  be  altogether  present  in  the  external  velocity  fields  calculated  by  the  constant 
source  approach  and  that  most  of  the  velocity  discontinuities  not  only  will  be  spurious  in  their 
own  right  but  also  will  produce  unwanted  repercussions  on  the  rest  of  the  configuration.  For 
such  method  straightforward  attempts  to  calculate  supersonic  flow  field  streamlines  are  likely  to 
produce  rather  odd  jagged  shapes  that  become  less  smooth  as  the  number  of  calculated  points 
increases. 

To  minimize  the  undesirable  effect  caused  by  the  internal  spurious  wave-train,  it  is  apparent  that 
techniques  to  improve  continuity  in  the  representation  of  both  geometric  and  surface  singularity 
variation  are  required.  The  PANAIR  Code  (Ref  3.30)  solves  the  steady  aerodynamics  by  a  high 
order  steady  panel  model.  Linear  source  and  quadratic  doublet  singularities  are  distributed  over 
each  body  box  and  the  continuity  conditions  of  singularity  distribution  across  adjacent  boxes  are 
imposed.  The  unknown  steady  source  strength  and  doublet  strength  are  solved  by  imposing 
Neumann  and  Dirichlet  boundary  conditions  on  each  boxes.  However,  extending  PANAIR’s 
high  order  panel  method  to  unsteady  one  is  extremely  complex  and  computationally  expensive. 

Instead  of  using  a  high  order  panel  method,  ZONA7  employs  the  so-called  “integral-averaged 
boundary  condition”  to  minimize  the  spurious  waves  resulted  from  the  constant  source 
singularity  (Ref  3.32).  The  integral-arranged  boundary  condition  is  a  wave  minimization 

technique.  It  defines  an  averaged  normal  velocity  Kn  on  each  box  by  a  sample  integral  average 
technique; 

F.  =  1  Jl  (V«(,  •  S)  (3.71) 

where  jp  denotes  the  box  area. 

For  steady  aerodynamics,  the  integral-averaged  boundary  condition  reads: 

F„=-(F-n)  (3.72) 
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Comparing  Eq  3.72  with  the  steady  boundary  condition  of  Eq  3.28,  it  can  be  seen  that  Eq  3.28  is 
normally  imposed  at  the  control  point  whereas  Eq  3.72  is  satisfied  over  the  box  area.  The 

averaged  normal  velocity  Fn  can  be  obtained  by  a  four-point  Gaussian  Quadrature  of  Eq  3.71. 
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Fig  3.40  Cone-Cylinder-Cone  at  Af  =  2.0  and  a  =  0°  Using  ZONA7  with  Wave  Minimization 

Fig  3.40  shows  the  results  obtained  by  ZONA7  with  the  integral-averaged  boundary  condition 
(wave  minimization  technique)  of  the  cone-cylinder-cone  case  at  M=  2.0.  Comparing  to  the 
one  without  wave  minimization  shown  in  Fig  3.37,  it  can  be  seen  that  the  stability  of  the  solution 
with  the  wave  minimization  technique  is  greatly  improved. 

It  is  found  that  the  integral-averaged  boundary  condition  works  equally  well  for  unsteady 
aerodynamics.  This  is  done  by  applying  Eq  3.71  to  the  unsteady  boundary  condition  of  Eq 
3.23b.  Thus,  the  integral-averaged  boundary  condition  offers  a  technique  to  minimizes  the 
spurious  wave  caused  by  the  constant  source  singularity  approach  for  both  steady  and  unsteady 
supersonic  aerodynamics. 

3.11  Super-Inclined  Boxes  and  Inlet  Boxes 

In  a  linearized  supersonic  flow  formulation,  the  freestream  Mach  cone  determines  the  region  of 
influence.  Typical  supersonic  panel  methods  generally  work  well  if  the  body  under 
consideration  is  fully  immersed  within  this  region  of  influence.  However,  when  the  supersonic 
freestream  becomes  higher  and/or  the  body  is  relatively  thick  whereby  a  part  of  the  body  would 
be  exposed  outside  of  the  zone  of  influence,  most  supersonic  panel  methods  would  cease  to  be 
applicable. 
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For  boxes  placed  on  the  inlet  surface  (Fig  3.41a)  or  on  the  nose  of  a  thick  body  (Fig  3.41b),  the 
local  angles  of  incidence  on  some  boxes  would  be  greater  than  the  freestream  Mach  cone  angle, 
this  would  render  them  lie  outside  of  the  freestream  Mach  cone.  These  boxes  are  called  “super- 
inclined  boxes”  and  they  are  the  causes  for  numerical  singularities  in  the  supersonic  aerodynamic 
influence  coefficient  computation. 


(a) 


(b) 

Fig  3.41  Superinclined  Box  (a)  on  Engine  Inlet  (b)  on  Thick  Body 

To  circumvent  this  numerical  singularity  problem  that  is  associated  with  super-inclined  boxes  in 
supersonic  flow,  we  introduce  a  special  treatment  for  the  aerodynamic  influence  coefficient 
computation  in  ZONA7.  This  engineering  treatment  adopts  the  corresponding  oblique  shock 
angle  for  a  cone  (based  on  the  Exact  Euler  Conical-Flow  Solutions,  e.g.,  Sims,  Ref  3.33, 3.34)  to 
substitute  the  Mach  wave  angle.  The  local  cone  angle  for  each  superinclined  box  is  measured  by 
the  angle  between  the  freestream  and  the  slope  of  the  box.  The  coiresponding  oblique  shock 
angle  is  used  as  a  modified  Mach  wave  angle  to  position  a  “Mach  Wave”  slightly  ahead  of  the 
super-inclined  boxes.  It  is  this  modified  Mach  angle  that  will  determine  the  region  of  influence 
of  the  super-inclined  boxes. 

3.12  Treatment  of  Engine  Inlets 

ZONA6/7  can  also  include  the  aerodynamic  computations  for  engine  inlets.  For  body  boxes  on 
the  inlet  face,  the  boundary  condition  must  be  modified  for  the  “flow-through”  condition.  This 
type  of  body  boxes  is  called  “inlet  box”. 
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Let  Ac  be  the  height  of  the  stream  tube  containing  the  maximum  mass  flow  which  can  enter  the 
inlet,  length  Ao  be  the  height  of  the  stream  tube  actually  entering  the  inlet,  the  mass  flow  ratio 
(MFR)  is  defined  as: 

MT?  =  (3.73) 

A 

For  the  inlet  operating  critically,  no  blockage  occurs  and  the  inlets  are  “flowing  full”.  In  this 
case  Ao  is  equal  to  Ac  and  AdFR  =  1.  For  the  subcritical  flow  condition,  only  the  stream  tube 
with  height  Ao  enters  the  inlet  The  remaining  stream  tube  with  height  (Ac  -  Ao)  is  spilled  from 
the  intake.  For  this  case,  MFR  becomes  less  than  1.  For  a  given  mass  flow  ratio  (MFR),  the 
steady  boundary  condition  expressed  in  Eq  3.28  can  be  modified  for  the  inlet  boxes,  as: 

V  .  K  =  -  (f  •  n)(l  -  MFR)  (3.74) 

and  the  unsteady  boundary  condition  becomes: 


=  Fb(1  -  AiF/?)  (3.75) 

Regions  OPEN  and  BLOCKED  are  determined  by  the  intersections  of  the  stream  tubes  and  the 
box  of  the  inlet.  As  such,  the  lengths  of  the  two  regions  are  proportional  to  the  heights  of  the 
stream  tubes,  Ao  and  (Ac  -  Ao),  respectively.  Thus,  the  ratio  of  the  inlet  area  covered  by  OPEN 
boxes  to  the  total  inlet  area  reflects  the  mass  flow  ratio  A^R,  where  MFR  is  determined  by  the 
intersection  of  the  stream  tubes  and  the  inlet  box. 

The  inlet  in  supersonic  freestream  usually  requires  the  treatment  of  superinclined  boxes  (see  Fig 
3.41a).  In  ZONA7,  the  calculation  of  Aerodynamic  Influence  Coefficients  (AIC)  for  an  inlet  box 
with  respect  to  other  boxes  adopts  the  modified  Mach  angle  substitution  developed  previously 
for  superinclined  boxes  of  thick  bodies. 
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4.0  ZTAIC:  UNSTEADY  TRANSONIC  AERODYNAMICS 
WITH  STEADY  PRESSURE  INPUT 


ZTAIC  is  ZONA’s  Transonic  Aerodynamic  Influence  Coefficient  method  for  generating 
unsteady  transonic  aerodynamics  for  lifting  surface  systems.  With  provided  steady  pressure 
distribution  from  CFD  solutions  or  measured  data  on  each  wing  section,  ZTAIC  generated 
transonic  (modal)  AIC’s,  which  necessarily  supports  effective  unsteady  aerodynamics/aeroelastic 
computations.  Hence,  in  addition  to  the  ZONA6  input,  the  supplied  vdng-section  steady 
pressures  are  required  to  be  inputted  by  the  user. 

Unsteady  subsonic  and  supersonic  aerodynamic  forces  and  moments  can  be  well  predicted  by 
ZONA6  and  ZONA7.  But  in  the  transonic  flight  regime  these  methods  cease  to  be  applicable 
due  to  the  nonlinear  compressibility  and  the  presence  of  transonic  shock  waves.  The  nonlinear 
Transonic  Small  Disturbance  Equation  (TSDE)  can  be  time-linearized  with  respect  to  the  motion 
amplitude.  This  allows  the  formulation  of  a  modal-based  ZONA’s  Transonic  Aerodynamic 
Influence  Coefficient  (TAIC)  matrix  method.  In  this  way,  ZTAIC  generates  unsteady  transonic 
aerodynamic  forces  and  moments  that  properly  account  for  the  transonic  effects  including  shock 
oscillations  up  to  the  linear  order  in  amplitude.  The  unique  features  of  ZTAIC  include: 

-  ZTAIC  input  requirement  is  nearly  identical  to  these  of  ZONA6  and  ZONA7  except  the 
additional  user-supplied  steady  pressure  input. 

-  The  user-supplied  steady  pressure  input  can  be  obtained  either  by  wind  tunnel  measured  data 
or  a  high  level  CFD  solver.  These  steady  input  contains  sufficient  detail  of  the  transonic 
shock  structure;  hence,  the  resulting  unsteady  pressures  in  terms  of  unsteady  shock  strength 
and  location  are  found  to  be  well  predicted. 

-  ZTAIC  does  not  require  grid  generation.  The  grid  system  for  solving  the  unsteady  transonic 
small  disturbance  equation  is  built  in  ZTAIC.  ZTAIC  can  automatically  optimize  the  mesh 
of  the  grid  system  according  to  given  Mach  number,  oscillating  frequency,  and  steady 
pressure  input. 

-  ZTAIC  generates  a  transonic  AIC  matrix  using  a  Modal-based  AIC  (MAIC)  approach. 
MAIC  represents  a  frequency-domain  aerodynamic  transfer  function  whose  definition  is  the 
same  as  that  of  ZONA6  and  ZONA7. 

In  this  section,  we  will  discuss: 

•  Background  of  ZTAIC 

•  Inverse  airfoil  design  scheme  based  on  the  users  supplied  steady  Cp  input 

•  Time-domain  unsteady  transonic  aerodynamic  computational  scheme  of  ZTAIC 

•  Transforming  ZTAIC  to  the  frequency-domain 

•  Modal  AIC  (MAIC)  approach 
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4.1  Background  of  ZTAIC 


In  the  transonic  flow  regimes,  the  wing  (body)  thickness  renders  the  problem  nonlinear  and  the 
mean  flow  nonuniform  (for  an  unsteady  problem).  Thus  the  current  treatment  of  the  unsteady 
transonic  mostly  adopts  the  Computational  Fluid  Dynamics  (CFD)  methodology.  Currently 
there  exists  a  number  of  well-practiced  CFD  unsteady  transonic  methods  presumably  ready  for 
aeroelastic  applications  (e.g.  CFL3D,  ENS3DAE  and  ENSAERO,  Refs  4.1,  4.2,  4.3).  However, 
their  acceptance  by  the  aerospace  industry  for  rapid  analysis  and  design  is  still  hampered  by 
problems  in  grid  generation,  CFD/CSD  interfacing  and  affordable  computing  time.  For  example, 
to  perform  a  high-level  (Navier-Stokes)  CFD  computation  for  a  single  wing  structure  flutter 
boundary  would  take  days  if  not  weeks  of  CPU  time  on  a  mainframe  computer.  Thus,  if 
rendering  such  high-level  CFD  methods  for  flutter  boundary  prediction  of  a  whole  aircraft 
should  take  up  much  longer  CPU  time  which  will  not  be  affordable  for  the  industry.  On  the 
other  hand,  solving  the  unsteady  TSDE  equation  using  a  lower-level  CFD  method  such  as  CAP- 
TSD  (Ref  4.4)  should  be  computationally  efficient.  But  CAP-TSD  is  confined  by  its  self¬ 
generated  steady-flow  solution  in  that  the  TSD  assumption  would  restrict  its  applicability  to 
wings  with  thick-airfoil  or  supercritical-wing  sections.  If  the  viscous  effect  is  not  introduced  in 
the  steady  mean  flow  of  an  unsteady  TSD  method  such  as  CAP-TSD,  then  the  resulting  shock 
position  and  strength,  and  hence  the  flutter  solutions  would  be  invalidated. 

Further  attempts  in  an  MDO  environment  have  been  somewhat  discouraged  by  their 
incompatibility  with  structural  FEM.  Although  current  CFD  efforts  are  directed  towards  meeting 
the  FEM  compatibility,  questions  of  affordable  computing  time  for  each  MDO  cycle  remain 
outstanding.  Toward  this  end,  we  have  re-examined  the  imsteady  transonic  aerodynamic 
methodology  critically  from  the  viewpoint  of  its  FEM  compatibility  and  its  expediency  for 
aeroelastic  analysis  and  design  applications.  The  result  of  this  re-examination  effort  is  the 
development  of  ZTAIC. 

Since  1988,  ZONA  has  been  following  up  on  the  development  of  the  Transonic  Equivalent  Strip 
(TES)  method  originally  supported  by  NAVAIR/ONR  (Refs  4.5,  4.6)  for  unsteady  flow 
computation  of  arbitrary  wing  planforms.  The  TES  method  consists  of  two  consecutive  steps 
added  to  a  given  nonlinear  Transonic  Small  Disturbance  code  such  as  ZTRAN  (Ref  4.7);  namely 
the  chordwise  mean  flow  correction  and  the  spanwise  phase  correction.  The  chordwise  mean 
flow  correction  is  accomplished  by  an  inverse  airfoil  design  scheme  incorporated  in  ZTRAN. 
The  computational  procedure  require  direct  steady  pressure  input  from  other  CFD  codes  or 
measured  data.  It  does  not  require  airfoil  shape  or  grid  generation.  Meanwhile,  all  the  mean- 
flow  shock  jumps  are  properly  included  in  the  resulting  unsteady  aerodynamics.  Based  on  the 
TES  concept,  ZONA  has  further  developed  the  ZONA  Transonic  Aerodynamic  Influence 
Coefficient  (ZTAIC)  code  by  adopting  a  Model-based  AJC  (MAIC)  procedure  ^ef  4.8).  ZTAIC 
has  been  successfully  applied  in  a  Multi-Disciplinary  Optimization  (MDO)  software  system  for 
aeroelastic  applications  (Ref  4.9).  Shown  in  Fig  4.1  is  a  flow  chart  presenting  ZTAIC  procedure 
for  flutter  analysis. 
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Fig  4.1  Flow  Chart  of  ZTAIC  Computation  Procedure 


Results  computed  by  ZTAIC  have  been  validated  with  wind  tunnel  data  for  a  number  of  wing 
planforms  with  various  aspect  ratios.  These  include; 

Unsteady  Pressures:  Measured  Data  Input 

•  The  Lessing  Wing  (Ref  4. 1 0) 

Fig  4.2a  shows  the  Lessing  Wing  with  a  rectangular  planform  of  aspect  ratio  3.0  and  width  a  5  % 
thick  parabolic  arc  airfoil  section.  Unsteady  pressure  magnitudes  and  phase  angles  are  presented 
in  Figs  4.2b  and  4.2c  for  spanwise  locations  at =  0.5,  0.7  and  0.9  (;;  =  2y  /  b).  It  is  seen  that 
the  ZTAIC  results  correlate  reasonably  well  with  two  sets  of  experimental  data  (1®*  and  2“*  runs). 

•  The  LANN  Wing  (Ref  4.11) 

Fig  4.3  shows  the  dimensions  of  the  LANN  Wing  Planform  of  aspect  ratio  7.93  and  a  12  %  thick 
airfoil  section  as  well  as  the  in-phase  and  out-of-phase  pressures  of  the  LANN  Wing  with 
pitching  axis  at  62  %  root  chord.  Throughout  three  spanwise  locations  considered  (7  =  32.5  %, 
65  %  and  95  %),  the  present  results  for  the  upper  surface  compare  more  favorably  with  NLR 
measured  data  that  do  the  XTRAN3  results  (Ref  4.12).  Since  subcritical  flows  are  predicted  for 
lower  surfaces,  the  unsteady  pressures  do  not  contain  the  shock  jump,  as  expected. 
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•  The  F-5  Wing  with  Control  Surface  (Ref  4. 13) 

Fig  4.4  compares  in  and  out-of-phase  pressures  at  two  flap  sections  of  F-5  wing  at  =  0.9. 
Close  agreements  are  found  with  XTRAN3S  (Ref  4. 14)  results  and  the  NLR  measured  data. 

Flutter  Results 

•  The  AGARD  Standard  445. 6  Wing  (Ref  4. 15) 

The  445.6  Wing  planform,  has  an  aspect  ratio  of  4  and  a  NACA  65 A004  airfoil  section.  It  has 
two  structural  models:  the  solid  wing  and  the  weakened  wing.  This  is  an  ideal  case  to 
demonstrate  ZTAIC’s  AIC  capability.  The  aerodynamic  shapes  of  these  two  models  remain  die 
same,  but  structurally  they  have  two  different  sets  of  structural  modes  associated  with  the  solid 
wing  and  the  weakened  vsdng.  For  this  reason,  the  same  modal  AIC  can  be  shared  by  both 
models.  Hence,  the  modal  AIC  computed  for  the  weakened  wing  can  be  saved  allowing  for  a 
warm-start  for  the  solid  iving.  Fig  4.5  presents  the  flutter  results  of  the  weakened  wing.  At  a 
subsonic  Mach  number  =  0.678,  the  ZTAIC  result  is  in  good  agreement  ivith  that  of 
ZONA6,  as  expected.  At  other  supercritical  Mach  numbers  where  =  0.9  and  0.95,  ZTAIC 
predicts  a  pronounced  transonic  dip  which  is  comparable  to  that  predicted  by  the  CAP-TSD 
code. 

Fig  4.6  presents  the  flutter  results  of  the  solid  wing.  The  stored  modal  AIC  from  the  weakened 
wing  is  warm-started  here.  Consequently  the  computing  time  for  each  flutter  point  requires  only 
one  minute,  in  contrast  to  the  weakened  wing  case  which  requires  five  hours  of  CPU  time  for  a 
flutter  point  on  a  SUN  SPARC  20  workstation. 
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Fig  4.5  Comparison  of  Flutter  Speed  and  Frequency  of  445.6  Weakened  Wing  at 
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Fig  4.6  Comparison  of  Flutter  Speed  and  Frequency  of  445.6  Solid  Wing  at  M  =  0.90  and  0.95 
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4.2  The  Chordwise  Mean  Flow  and  Spanwise  Phase  Correction  Procedures 

Governin2  Equations 

The  unsteady  Transonic  Small  Disturbance  (TSD)  equation  reads: 


[(1  -  Ml)  -  (r  +  l)Ml  +  4>„  +  «>„  -  =  0  (4.1) 

Similar  to  Eq  3.2,  the  TSD  potential  <J>  can  be  expressed  in  terms  of  the  steady  TSD  potential 
and  the  time-linearized  TSD  potential  with  respect  to  linear  amplitude,  i.e.: 

^  +  ^{x,y,z,t)  (4.2) 

where  and  should  satisfy  the  following  equations,  respectively. 

[(1  -  Ml)  -  (r  +  l)Ml  =  0  (4.3) 

(1  -  +  (y  +  +  ^yy  +  =  2Ml^„  +  Ml^„  (4.4) 

Eq  4.3  is  the  famous  steady  TSD  equation  due  to  von  Karman.  Eq  4.4  is  the  so-called  time- 
linearized  equation  derived  by  Landahl  (Unsteady  Transonic  Flow,  Pergamon  Press,  1961,  Ref 
4.16). 

Notice  that  Eqs  4.3  and  4.4  are  the  transonic  counterparts  of  Eqs  3.4  and  3.5  for 
subsonic/supersonic  linearized  formulation. 

The  premises  of  the  Transonic  Equivalent  Strip  (TSE)  method  (Ref  4.5)  is  to  propose  utilizing  a 
two-dimensional  TSD  potential  solver  to  solve  Eq  4.4  in  conjunction  with  a  3D  linear  potential 
solver  such  as  ZONA6.  LTRAN2  (Ref  4.17)  was  originally  selected  for  this  2D  solver.  But 
ZONA'S  improved  ZTRAN  is  to  replace  LTRAN2  in  all  our  subsequent  development.  The 
governing  2D  TSD  equation  is  nothing  more  than  a  2D  version  of  Eq  4.1,  viz.: 

Ui-Ml)  -  (r+l)W>J®„  +  (4.5) 

The  LTRAN2  of  Ballhaus  and  Gooijian  solves  Eq  4.5  with  two  alternatives:  the  Harmonic 
Oscillation  Method  and  the  Indicial  Method  (Refs  4.17, 4.18).  The  latter  method  amounts  to  an 
expedient  method  of  time  linearization  of  Eq  4.5  directly  by  CFD  means.  Equivalently,  their 
indicial  method  solves  a  2D  version  of  Eq  4.4,  viz.: 

(l~Ml)(p^^  +  ir  +  =  2Ml<p,^  +  (4.6) 

where  ,  the  2D  steady  TSD  potential,  satisfies  the  2D  Karman  equation,  i.e.: 
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-  (r +  !)«>„  ]?>™,  +  «>„  =  0 


(4.7) 


Boundary  Conditions 

On  the  mean  surface  of  the  wing  planform,  the  potentials  must  satisfy  the  linearized  tangency 
condition,  i.e.  at  z  =  0, 


rT{x,y„0) 

(4.8) 

S[H^  +  H,} 

where:  T  is  the  wing  airfoil-section  thickness  distribution  at  spanwise  station  y\ 

H  is  the  mode  shape 

r  is  the  thickness  ratio 

5  is  the  oscillation  amplitude,  and  6  <  r. 

The  tangency  condition  along  with  the  wake  condition  imposed  here  are  the  typical  linearized 
boimdary  conditions,  also  adopted  in  the  previous  linearized  subsonic  and  supersonic  potential 
formulations.  Hence,  one  assumes  that  4>o  ^re  decoupled  throughout  except  the 

governing  equation,  Eq  4.4,  itself.  For  TSD  potentials  <p^  and  of  Eqs  4.6  and  4.7,  they 
should  satisfy  the  same  boundary  conditions  as  above  except  in  the  2D  sense. 

TES  Solution  of Ea  4. 4 

Proposed  earlier  by  Liu  et  al.  (Ref  4.6)  and  further  supported  by  Oyibo’s  Separability  Principle 
for  Full-Potential  Transonic  Solutions  (Ref  4. 19),  a  solution  of  Eq  4.4  can  be  derived,  i.e.,  at 

a  spanwise  location  y  =  yu 

^i(x,yi,z,/)  =  (Z),,(x,z,0  F(yi;Ai)  (4.9) 

where:  is  the  2D  unsteady  solution  by  ZTRAN  at  a  spanwise  location  y  =  yi 

F  is  Oyibo’s  spanwise  exponential  decaying  function 

Ai  is  the  decaying  parameter  measured  by  the  chordwise  (2D)  to  the  spanwise 
(3D)  lift  slope  ratio. 

Clearly,  the  same  separability  of  in  Eq  4.9  showed  also  holds  for  of  Eq  3.5,  i.e.: 

^/(jr,yi,z,0  =  (py(x,z,t)  G(yi;«i)  (4.10) 

where  the  superscript  “  e  ”  denotes  the  purely  “linear”  subsonic  perturbed  potential. 

Pressure  coefficients 

The  linearized  unsteady  pressure  coefficient  reads: 


73 


AC,  =  +4) 


(4.11) 


Again,  Eq  4.1 1  holds  in  general  for  unsteady  transonic  potential  or  q>^  as  well  as  for  unsteady 

subsonic  potential  or  Combining  Eqs  4.9,  4.10  and  4.11  yields  the  TES  pressure 
relations  for  the  correction  procedure,  i.e.  at^  =  y\. 

AC,/  =  AC,/  /„(  AC,/,AC,/ )  .  exp{(;ii  }  (4.12) 

where  superscripts  N  and  t  denote  transonic  and  subsonic  unsteady  pressures  and  AC,/  and 

AC,/  denote  the  3D  and  2D  imsteady  pressure,  respectively.  Hence,  the  proposed  correction 
procedures  for  TES  lies  in  the  pressure  relation  given  by  Eq  4.12.  The  chordwise  mean-flow 
correction  is  supplied  by  “injecting”  the  correct  steady  pressure  input  in  AC,/ ,  which  is  to  be 

computed  by  ZTRAN.  Recall  that  steady  pressure  input  is  provided  by  measured  data  or  by  a 
high-level  CFD  solution.  The  way  to  “inject”  this  input  into  ZTRAN  is  accomplished  by  an 
inverse  airfoil  design  (IAF2)  scheme  to  be  discussed  in  the  next  section.  The  spanwise  phase 

correction  is  accomplished  by  providing  the  linear  pressure  distributions  AC,^'  and  in 

the  function  “/» “,  whereas  they  are  obtained  from  ZONA6  computation. 

The  exponential  decaying  parameters  A,  and  are  measures  of  3D  to  2D  lift  curve  slope  ratios 
(for  transonic  and  subsonic,  respectively).  For  simplicity,  the  current  TES  procedure  maintains 
at  an  approximated  level  where  A,  =  a^. 

Finally,  the  two  correction  procedures  combined  in  Eq  4.12  have  the  physical  meaning  that:  i) 
the  chordwise  correction  accounts  for  reproducing  the  nonlinear  structure  of  the  three- 
dimensional  mean  flow,  ii)  the  spanwise  phase  correction  is  responsible  for  the  adjustment  of 
the  spanwise  phase  lag  of  the  pressure  according  to  an  equivalent  linear  three-dimensional  flow. 
Clearly,  shock  waves  caimot  be  created  or  destroyed  by  any  process  of  these  corrections. 

4.3  The  Inverse  Airfoil  Design  (1AF2)  Scheme 

It  has  been  pointed  out  by  Fung  (Ref  4.20)  and  Lamboume  (Ref  4.21),  among  others,  that  an 
accurate  steady  state  vith  correct  shock  jump  and  locations  is  essential  for  correct  unsteady 
aerodynamic  computations.  It  is  believed  that  TES  in  general  should  be  adequate  in  handling  the 
classical  flutter  analysis  in  the  transonic  regime,  because  its  chordwise  mean-flow  correction 
could  result  in  the  proper  shock  position  and  strength  including  viscous  effects.  This  correction 
procedure  is  made  possible  by  adopting  the  externally  provided  steady  pressure  input  from 
measured  data  or  from  high-level  CFD  computations.  To  “inject”  the  steady  pressure  input  into 
ZTRAN  of  the  TES  method  requires  to  build  in  an  equivalent  airfoil  design  procedure  using  an 
inverse  design  algorithm  called  IAF2.  In  doing  so,  a  proper  steady  flowfield  can  be  generated 
corresponding  to  the  input  steady  pressure  prescribed  on  the  airfoil  surface. 
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(A)  lANN  WIilC 
M_  -  0.82 


(B)  HOKTHROf  F-S  WINC 
M_  -  0.90 


(C)  HORTHSOP  F-S  WIHC 
M_  -  0.95 


Fig  4.7  Steady  Pressure  Inputs  and  Equivalent  Airfoil  Outputs  at  various  Spanwise  Locations 

(□  upper  surface,  A  lower  surface, - presents  TES  data): 

a)  LANN  Wing  at  Mean  Incidence  oc<,  =  0.62®,  Af„  =  0.82  (NLR  Measured  Data); 

b)  Northrop  F-5  Wing  at  Mean  Incidence  oc^  =  0°,  =  0.90  (NLR  Measured  Data); 

c)  Northrop  F-5  Wing  at  Mean  Incidence  oco  =  0®,  =  0.95  (XTRANSS-Ames  Input) 
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Based  on  the  user  supplied  steady  pressure  input,  an  Inverse  Approximate  Factorization  scheme 
(IAF2)  proposed  by  Fung  and  Chung  (Ref  4.7)  is  adopted  to  perform  an  inverse  airfoil  design 
according  to  Eq  4.7.  The  velocity  potential  obtained  from  integrating  the  pressure  on  the  slit 
representing  an  airfoil  is  known  up  to  an  arbitrary  constant.  To  determine  this  constant,  it  is 
required  that  a  closure  condition  is  imposed,  e.g.  the  resulting  slope  distribution  being  equivalent 
to  a  closed  body.  This  constant  a  being  updating  during  the  numerical  iteration  process  until  a 
converged  solution  satisfying  the  closure  requirement  is  obtained.  Fig  4.7a  and  Fig  4.7b  display 
the  pressure  inputs  for  upper  and  lower  surface  at  0.82  and  0.9  for  the  LANN  and  F-5 
wings,  respectively,  based  on  the  NLR  measured  data.  The  input  data  in  general  are  shown  by 
open  symbols  whereas  the  pressures  computed  by  the  IAF2  scheme  by  solid  lines. 

Fig  4.7c  demonstrates  the  flexibility  of  the  pressure  input  scheme  using  inverse  airfoil  design,  in 
which  the  computed  pressures  by  XTRAN3S  of  F-5  wing  at  M=0.95  is  adopted  as  the  steady 
pressure  input.  The  strong  shock  close  to  the  trailing  edge  in  this  case  is  considered  a 
challenging  case  for  an  airfoil  design  code.  As  can  be  seen  in  Fig  4.7c,  the  present  IAF2  showed 
excellent  comparison  between  the  XTRAN3S  pressure  input  and  the  design  output.  Further 
studies  in  using  the  present  IAF2  for  pressure  input/output  verification  showed  the  present 
scheme  is  indeed  very  robust  and  can  generate  accurate  pressure  flowfield  associated  with  a  wide 
class  of  transonic  pressure  input  with  transonic  shock  jumps. 

4.4  Frequency-Domain  Pressure  Coefficient  by  Indicial  Method 

For  2D  unsteady  transonic  computation,  the  ZTRAN  code,  similar  to  LTRAN2,  also  has  two 
alternative  methods,  namely  the  harmonic  oscillation  method  and  a  method  based  on  a  time 
linearized  scheme  proposed  by  Fung  et  al.  (Ref  4.22).  Most  unsteady  TSD  codes  like  LTRAN2 
and  XTRAN3S  generated  their  steady  aerodynamics  and  pressures  according  to  the  steady 
TSDE.  Unlike  these  methods,  ZTRAN  employs  the  inverse  airfoil  design  to  inject  the  externally 
provided  steady  pressure  input  and  generate  Ae  corresponding  steady  computation.  Thus,  this 
steady  pressure  input  procedure  removes  the  mandatory  usage  of  self-generated  steady 
aerodynamics.  Instead  the  steady  pressure  input  can  be  supplied  from  the  results  of  a  high-level 
CFD  computation  or  wind  tunnel  measured  data.  In  this  way,  the  inverse  airfoil  design  scheme 
ensures  that  the  nonlinear  effects  including  shock  structure  and  viscous  effects  are  closely 
reproduced  and  properly  incorporated  into  the  follow-on  unsteady  aerodynamic  computation 
scheme. 

Both  the  harmonic  oscillation  scheme  and  the  time-linearized  scheme  of  ZTRAN  yield  the 
unsteady  aerodynamic  solution  in  time  domain.  The  former  scheme  can  be  further  linearized  if 
the  oscillatory  amplitude  is  kept  sufficiently  small. 

To  achieve  a  sinusoidal  oscillation  solution  (linearized),  it  is  required  to  apply  the  following 
procedures.  First,  the  aerodynamic-response  solution  is  to  be  selected  only  after  a  few  cycles  of 
oscillation,  when  the  transient  solution  is  seen  to  be  diminishing  (or  the  response  solution 
becomes  closely  periodic).  Next,  Fourier  analysis  is  applied  to  the  selected  last-cycle  response 
solution  in  which  only  the  first  harmonic  component  is  retained,  whereas  all  the  higher 
harmonics  are  filtered  through  in  this  step.  From  the  first-harmonic  solution,  the  pressure 
magnitude  and  phase  angle  can  be  properly  identified.  This  method  for  the  first-harmonic 
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solution  works  fine  for  the  K  and  P-K  flutter  methods,  except  that  each  solution  obtained 
corresponds  to  only  one  single  frequency  input,  hence  it  tends  to  be  computationally  ineffective. 
For  computational  expediency,  the  indicial  method  is  preferred. 

The  indicial  method,  first  introduced  by  Tobak  (Ref  4.23),  offers  a  rapid  procedure  for 
converting  the  time-domain  solution  of  a  CFD  code  to  the  frequency  domain.  Normally,  the 
objective  of  the  indicial  method  is  to  obtain  the  “aeroelastic  transfer  function”  of  an  aerodynamic 
quantity  of  interest  in  the  frequency  domain,  i.e. : 

B.iik)  =  (4.13) 

where: 

-ikT 

?(.)  =  £  (Oe  dT  (4.14) 


is  the  Fourier  transform  of  (•),  T  =  is  the  non-dimensional  time  and  H  denotes  the 

L 

aeroelastic  transfer  function  of  an  aerodynamic  quantity  of  interest  such  as  Cp,  Cl,  ...  etc. 


Ballhaus  and  Gooijian  (Ref  4.24)  have  demonstrated  the  applicability  of  the  indicial  method  for 
transonic  aerodynamics.  However,  the  step  function  of  the  indicial  method  has  numerical 
difficulties  due  to  the  evaluation  of  the  Dirac  delta  function  associated  with  the  derivative  of  the 
step  function  at  r  =  0.  Later,  Seidel,  Benett  and  Ricketts  (Ref  4.25)  introduced  a  pulse  function 
method  in  the  XTRAN3  code.  The  pulse  function  is  defined  as: 


f{T)  =  (4.15) 

where  Tc  is  the  non-dimensional  time  at  which  the  maximum  amplitude  is  reached  and  W 
represents  the  width  of  the  pulse.  The  desired  aeroelastic  transfer  fimction  is  then  obtained  by 
solving  flie  following  equation  numerically: 


H(/it) 


9  mi)) 

9  m)) 


(4.16) 


where  X(T)  is  the  aeroelastic  response  subject  to  the  pulse  function^?). 

Because  the  pulse  function  and  its  derivative  are  continuous,  it  removes  the  numerical  difficulty. 
Eq  4.16  also  suggests  that  the  entire  frequency-domain  solution  can  be  obtained  by  Fourier 
Transform  on  X{T)  and  f(J)  in  a  post-processing  sense.  Therefore,  instead  of  repeated 
computation  for  each  frequency  by  the  sinusoidal  method,  the  pulse  function  method  provides 
the  entire  frequency-domain  unsteady  aerodynamic  solution  with  only  a  single  imsteady 
aerodynamic  computation,  leading  to  a  rapid  procedure  for  frequency-domain  solution. 
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4.5  Modal  Aerodynamic  Influence  Coefficient  (MAIC)  Approach 

A  formal  AIC  matrix  should  contain  purely  aerodynamic  information  relevant  to  the  governing 
equation.  However,  to  generate  this  type  of  AIC  from  CFD  methods  such  as  the  ZTRAN  code 
remains  a  major  undertaking.  The  validity  of  the  MAIC  approach  is  based  on  the  assumption  of 
the  amplitude  linearization  principle  that  states:  “//ze  linearization  of  the  aerodynamics  for  an 
aeroelastic  system  in  any  flow  regime  can  be  assured  if  the  modal  amplitude  is  kept  stifficiently 
small  at  all  times'*.  In  fact,  this  is  equivalent  to  the  principle  of  time  linearization.  This  principle 
suggests  that  an  expedient  modal-based  AIC  procedure  can  be  developed  readily. 

Consider  a  typical  CFD  computation  procedure,  in  which  the  incremental  pressure  distribution  is 
related  to  a  given  pre-defined  baseline  modal  vector  by: 

{Cp},  =  /=l,...,m  (4.17) 

where  A^is  a  nonlinear  CHD  operator  that  generates  the  incremental  pressure  due  to  a  given 
baseline  modal  vector  (f) .  For  m  number  of  baseline  modal  vectors,  Eq  4.17  becomes  a  matrix 
equation: 

[P]  =  A'EH]  (4.18) 

where: 

m  =  [{C,},  ...  {C,},]  (4.19) 

is  defined  as  the  incremental  pressure  matrix. 

and: 

m  =  m,  (4.20) 

is  defined  as  the  baseline  modal  matrix. 

With  the  amplitude  linearization  principle  imposed,  the  operator  A' can  be  approximated  by  a 
linearized  operator  Thus: 

[P]  =  Xm  (4.21) 

The  baseline  modal  matrix  is  defined  as  a  set  of  displacement  vectors  that  could  be  superimposed 
to  closely  represent  the  deformation  of  the  structure.  Let  h  be  the  structural  deformation  vector, 
the  relation  between  h  and  the  baseline  modal  matrix  reads: 

{h}  =  [a>]{a}  (4.22) 

where  a’s  are  the  best-fit  coefficients,  to  be  determined  by  the  following  least  squares  procedure: 
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{a)  =  [([®f[4>]]’[«>]MW  ('*■23) 

Let  Cp  be  the  incremental  pressure  vector  due  to  the  structural  deformation  h,  then  Cp  can  be 
related  to  h  by  the  same  linear  operator: 

{C,)  =  Z{h}  (4.24) 

Substituting  Eq  4.22  into  Eq  4.24  yields; 

{Cp}  =  xm{a}  (4.25) 

Combining  Eqs  4.21,  4.23  and  4.25  gives  a  matrix  that  directly  related  the  incremental  pressure 
vector  to  the  structural  deformation  h: 

{Cp}  =  [MAlC\{\i)  (4.26) 

where: 

[MA1C\  =  [  [P]  [[Of  [0]]-‘  [O]^  ]  (4.27) 

is  the  sought  modal-based  44/C  matrix. 


The  modal-based  AlC  approach  is  a  general  procedure  which  could  be  applied  to  any  CFD  codes. 
Chen  et  al  demonstrated  the  MAIC  approach  in  a  transonic  aeroelastic  optimization  problem 
(Ref  4.9).  In  Ref  4.9,  the  structural  modes  of  the  initial  structure  are  defined  as  the  baseline 
modal  matrix.  The  modal-based  AIC  matrix  is  then  constructed  from  the  baseline  modal  matrix 
and  their  associated  unsteady  pressure  matrix.  The  constructed  modal-based  matrix  is  computed 
once  and  reused  in  the  structural  optimization  loop.  In  so  doing,  considerable  saving  in 
computer  time  is  achieved. 

The  ZTAIC  method  requires  a  somewhat  different  type  of  MAIC  matrix.  For  a  stripwise 
unsteady  transonic  aerodynamic  computation,  only  a  2D-baseline  modal  AIC  matrix  is  needed. 
For  high  aspect  ratio  wing  structures  involving  primarily  spanwise  bending  and  torsion 
deformations,  2D  rigid  pitch  mode  and  2D  rigid  plunge  mode  are  used  to  represent  the  local 
structural  deformation  at  each  strip.  For  wing  sections  containing  leading  edge  and  trailing  edge 
control  surfaces,  a  2D  leading  edge  and  2D  trailing  edge  modes  are  used  to  represent  these 
control  surfaces  at  the  strip.  For  low  aspect  ratio  wing  planforms  where  chordwise  bending 
deformation  may  occur,  a  2D  chordwise  bending  mode  is  introduced.  Altogether,  five  modal 
vectors  are  used  to  define  the  2D  baseline  modal  matrix. 

With  the  2D  baseline  modal  matrix  in  hand,  the  generation  of  the  MAIC  matrix  for  ZTAIC  is  a 
straight  forward  procedure.  The  formal  AIC  for  relating  the  structural  deformation  to  the 
aerodynamic  forces  can  be  obtained  simply  by  multiplying  the  area  of  each  aerodynamic  boxes 
to  the  MAIC  matrix.  This  gives: 
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{F,}  =  qAAIC\{h} 


(4.28) 


where: 

[AIC\  =  [n\[AREA]{MAIC] 


and: 

{AREA[\  is  a  diagonal  matrix  containing  the  area  of  each  aerodynamic  boxes. 

[«]  is  a  K-set  by  J-set  matrix  containing  the  normal  vector  components  of  each 
aerodynamic  box.  The  exact  expression  of  [«]  can  be  found  in  Eq  3.62. 
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5.0  ZONA7U:  UNIFIED  HYPERSONIC/SUPERSONIC 
UNSTEADY  AERODINAMICS  FOR  WING-BODY 
AIRCRAFT  CONFIGURATION 


ZONA7U  extends  the  applicability  of  the  lifting  surface  method  of  ZONA7  into  the  unified 
hypersonic  and  supersonic  flight  regime.  The  unsteady  lifting  surface  option  of  ZONA7U 
(formerly  ZONA51U)  takes  into  account  the  nonlinear  wing  thickness  effect  and  the  flow 
rotationality  effect  due  to  strong  shock  waves,  whereas  these  effects  are  neglected  by  the  linear 
theory  and  overestimated  by  Piston  theory.  In  addition  to  the  ZONA7  input,  only  the  wing- 
section  profiles  are  required  to  be  inputted  by  the  user.  In  this  section,  we  will  discuss: 

•  Background  of  ZONA  7  U 

•  Review  of  Piston  theory  and  Third-Order  theories 

•  Hypersonic  Similitude 

•  Unified  Hypersonic/Supersonic  Lifting  Surface  method  of  ZONA  7 U 

•  AIC  matrix  of  ZONA7U 

5.1  Background  of  ZONA7U 

Exact  unsteady  supersonic  three-dimensional  (3D)  linear  theory  has  been  sufficiently  developed 
for  the  treatments  of  lifting  surfaces,  e.g.  the  ZONA7  method.  Nonetheless,  the  lifting  smface 
method  incorporated  in  ZONA7  is  confined  to  flat-plate  wing  sections,  in  which  no  thickness 
effect  can  be  accounted  for.  But  such  a  thickness  effect  at  moderate  to  hi^  supersonic  Mach 
number  is  of  practical  importance,  for  it  should  usually  lower  the  predicted  flutter  speed. 

The  lifting  surface  method  of  ZONA7U  adopts  the  imified  supersonic  and  hypersonic  lifting 
surface  method  ZONA51U  (Ref  5.1).  ZONA51U  combines  the  supersonic  lifting  surface 
method  of  ZONA51  with  a  nonlinear  correction  matrix  Ey  based  on  Donov  &  LiimelTs 
imiformly-valid  higher-order  hypersonic/supersonic  scheme.  For  aeroelastic  applications, 
ZONA51U  has  been  applied  to  various  wing  planforms  with  thickness  distributions.  Some  of 
these  are  shown  as  follows: 

•  Rectangular  Wing  with  Thickness  Profile 

Fig  5.1  presents  the  variation  of  generalized  aerodynamic  forces  (C^^ )  with  reduced  frequency  k 

for  an  oscillating  wedge  ( cr  =10  )  at  M=  3.0.  Results  of  ZONA5 1  and  ZONA5  lU  are  compared 
with  that  of  the  Euler-Perturbed  Euler  Characteristics  (PEC)  method  (Ref  5.2).  Good  correlation 
is  found  between  the  results  of  ZONA51U  and  the  Euler-PEC  solution.  Substantial  departures 
exist  between  the  ZONA51  and  ZONA51U  results  throughout  the  frequency  range,  indicating 
clearly  the  persistent  nonlinear  thickness  effect. 
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Fig  5.2  Damping-in-pitch  vs.  Semiwedge  Fig  5.3  Stiffness  Derivative  C  „  ^  and  Damping-in- 

Angle,  A  =  0.5c:  M=  a)  5.0,  b)10.0  pitch  Derivatives  vs.  Semiwedge  Angle, 

h  =  0.5c:  M=  a)  5.0,  b)  10.0 

Fig  5.2  presents  the  variation  of  damping-in-pitch  derivative  with  wing  thickness  cr  for  a 
diamond  airfoil  section  at  2.0,  5.0  and  10.0.  For  M=  2.0,  nonlinear  results  of  Van  Dyke 
(Ref  5.3),  Piston  theory  and  ZONA51U  agree  well  up  to  cr  =15°.  When  Mach  number  is 
increased  to  M=  5.0  and  10.0,  all  nonlinear  results  tend  to  over-predict  the  damping  derivatives 
except  that  of  ZONA5 lU.  In  fact,  5  lU  yields  results  in  close  agreement  with  Hui’s  Exact  Euler 
Solution  (Euler-Hui,  Ref  5.4)  up  to  cr  =15°  for  all  Mach  numbers  considered.  By  contrast,  the 
linear  theory  such  as  ZONA51  fails  to  account  for  the  thickness  effect. 
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Fig  5.3  presents  the  Mach  number  variation  of  the  stiffness  derivative  and  damping-in-pitch 
derivative  for  a  diamond  profile  with  thickness  (a  =15").  It  is  seen  that  the  ZONA51U  results 
essentially  follow  the  trend  of  the  Euler-Hui  solutions  throughout  the  Mach  range.  Other 
nonlinear  results  begin  to  depart  from  the  Euler  solution  around  M=  4.0  and  continue  to  diverge 
from  it  as  the  Mach  number  increases.  In  fact,  among  all  of  the  above  results,  only  results  of 
ZONA51U  and  the  Euler-Hui  exact  solution  would  correctly  approach  the  Newtonian  limit. 


•  Panel  Flutter  ; 

Shown  in  Fig  5.4  are  two  flexible  panels  (membranes)  mounted  on  both  surfaces  of  a  wedge 
( or  =2").  The  motion  of  the  oscillating  panels  is  described  by  the  assumed  mode: 


e 


iki 


where  A"  =  2  and  s  is  the  amplitude  of  vibration. 


5.4  Oscillating  Panels  Mounted  on  a  Wedge  with  Semi-Wedge  Angle  a  =  2° 

Fig  5.5  presents  the  effect  of  reduced  frequency  on  generalized  aerodynamic  forces  for  these 
vibrating  panels  at  M=  5.0.  Generalized  Aerodynamic  Forces  (GAF)  results  in  Qn,  based  on  the 
linear  theory.  Piston  theory  and  ZONA51U  are  compared  with  the  Euler-PEC  solution.  Sirnilar 
to  the  earlier  observation  in  the  case  of  flap  oscillation,  ZONA51U  in  Fig  5.5  substantially 
improves  the  pressure  magnitude  over  that  of  the  linear  theory.  For  the  reason  stated  earlier, 
only  slight  improvement  of  the  phase  change  is  found  in  the  present  result.  Once  again,  the 
results  of  Piston  Theory  show  practically  no  phase  change  for  all  cases  considered. 
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Fig  5.5  Effect  of  Reduced  Frequency  on  Generalized  Aerodynamic  Forces  for  Oscillating 

Panels  {M  =  5.0,  a  =  2®,  iV=  2) 


•  SAAB  Canard-Wing 


Shown  in  Fig  5.6  is  the  aerodynamic  model  of  a  SAAB  Canard-Wing  combination.  Biconvex 
airfoil  sections  of  6%  in  thickness  are  assumed  for  both  canard  and  wing.  Two  modes,  the  wing 
plunging  and  the  canard  pitching  about  its  midchord,  are  considered. 


Fig  5.6  Paneling  Scheme  of  SAAB  Canard-Wing 
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In  Fig  5.7,  variation  of  generalized 
aerodynamic  forces  Qn  with  Mach 
number  is  presented  at  three  reduced 
frequencies,  k  =  0.01,  0.5  and  2.0.  It 
is  seen  that  the  trends  of  Qn  due  to 
ZONA51U  is  different  from  that  of 
the  linear  theory  in  the  hypersonic 
range  (between  M  =  5.0  and  M  = 
20.0)  showing  strong  canard-wing 
aerodynamic  interaction.  This  range 
is  found  to  be  related  to  the  canard- 
wing  dimensions.  As  can  be  seen 
from  Fig  5.6,  such  interaction  is 
expected  to  vanish  when  the 
hypersonic  Mach  number  is  increased 
up  to  M  =  21,  where  the  first  Mach 
wave  of  canard  detaches  from  the 
wing  trailing  edges. 
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Fig  5.7  GAF  Qu  of  SAAB  Canard-Wing  vs.  Mach  Number: 
Mode  1  Wing  Plunging  and  Mode  2  Canard  Pitching 
about  Midchord 


•  Wing  Flutter 

Two  wing  planforms  are  selected  for  performing  flutter  analysis  using  ZONA51U: 
a  70-degree  delta  wing  and  a  15-degree  swept  untapered  wing. 

•  70-Degree  Delta  Wing 

Fig  5.8a  presents  flutter  boundaries  for  a  70-degree  delta  wing  with  a  6%  thick  diamond  airfoil 
section  computed  by  ZONA5  lU  and  Piston  theory. 

The  flutter  experiment  was  carried  out  at  NASA  Langley  by  Hanson  and  Levey  (Ref  5.5).  The 
wing  model  used  was  essentially  a  flat-plate.  According  to  Ref  5.5,  four  measured  modes  are 
used  in  the  present  flutter  analysis.  Half  of  the  delta  planform  is  subdivided  into  10x10  panels. 
The  flutter  boundary  consists  of  the  flutter  points  obtained  for  six  Mach  numbers  (M=  1.19, 
1.30,  1.64,  2.0,  2.25  and  3.0).  Flutter  results  computed  by  CAP-TSD  (Ref  5.6),  Piston  theory, 
ZONA51  and  ZONA51U  are  compared  with  the  measured  data.  Several  observations  on  Fig 
5.8a  can  be  put  forth; 

-  ZONA51U  predicts  flutter  boundary  that  is  more  conservative  than  that  of  ZONA51 
indicating  that  the  thickness  effect  indeed  reduces  the  supersonic  flutter  speed.  Piston 
theory  on  the  other  hand  tends  to  over-predict  the  flutter  speeds. 

-  In  fact,  ZONA51U  predicts  the  most  conservative  flutter  boundary  among  all  methods 
considered.  Flutter  boundary  of  CAP-TSD,  on  the  other  hand,  appears  to  be  the  most 
non-conservative. 
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Fig  5.8a  Flutter  Speeds  and  Frequencies  vs.  Mach  Number  Predicted  by  Various  Mefliods:  70 

Delta  Wing  widi  and  without  Thickness 


Fig  5.8b  Flutter  Speed  and  Flutter  Frequency  for  a  15°  Swept  Untapered  Wing 

(M=  1.3  and  3.0) 
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15-Degree  Swept  Untapered  Wing 


Fig  5.8b  presents  two  computed  flutter  points  for  a  15-degree  swept  imtapered  wing  of  aspect 
ratio  AR  =  5.35  at  M=  1.3  and  M-  3.0. 

The  flutter  experiment  was  carried  out  at  NASA  Langley  Field  by  Tuovila  and  McCarty  (Ref 
5.7).  The  wing  model  used  is  a  cantilever  wing  with  a  2%  thick  hexagonal  airfoil  section. 
According  to  Ref  5.7,  eight  modes  generated  by  MSC/NASTRAN  are  used  in  the  present  flutter 
analysis.  Half  of  the  wing  planform  is  subdivided  into  10x10  panels. 

Computed  results  of  ZONA51,  Rodden’s  method  (Ref  5.8)  (employing  ZONA51),  and 
ZONA51U  are  compared  with  test  data  of  Tuovila  and  McCarty.  Note  that  Rodden’s  method 
adopts  coefficients  from  Van  Dyke’s  theory  and  Piston  theory  and  corrects  upon  results  of 
ZONA51.  While  the  predicted  flutter  speeds  due  to  Rodden  and  ZONA51U  are  slightly  non¬ 
conservative  at  M=  1.5,  both  are  conservative  at  M=  3.0.  Overall,  the  listed  results  confirms 
once  again  the  impact  of  thickness  on  flutter  speed.  The  linear  theory,  as  computed  by  ZONA51, 
yields  non-conservative  flutter  points  at  both  Mach  numbers. 

5.2  Review  of  Piston  Theory 

Subsequent  to  the  original  publication  of  Lighthill  (Ref  5.9),  Ashley  and  Zartarian  (Ref  5.10) 
first  proposed  the  application  of  Piston  Theory  for  flutter  analysis  and  other  aeroelastic 
applications.  They  found  that  the  nonlinear  thickness  effect  provided  by  the  theory  indeed 
results  in  a  more  conservative  flutter  boundary,  which  was  validated  by  measured  data.  Based 
on  a  criterion  that  if  any  one  of  the  conditions  holds,  namely: 

M^»l,  kM^»\  or  (5.1) 

Landahl,  Ashley  and  Mello-Christiansen  (Ref  5.11)  further  established  a  consistent  linearized 
Piston  Theory.  With  this  theory,  they  obtained  an  explicit  flutter  solution  for  a  typical  two 
dimensional  wing  section.  The  flutter  speed  according  to  their  theory  approaches  those  predicted 
by  the  exact  linear  theory  (Ref  5.12)  as  the  Mach  number  increases,  whereas  they  tend  to  depart 
from  the  latter  as  the  Mach  number  decreases  toward  unity. 

Originally,  Lighthill’s  Piston  Theory  accounts  for  the  effect  of  the  nonlinear  thickness  in  the  high 
Mach  number  range  such  that  » 1 .  It  imposes  the  condition  that  the  magnitude  of  the 
piston  velocity  never  exceeds  the  speed  of  sound  in  the  undisturbed  fluid.  The  aerodynamics  of 
this  analogy  is  to  state  that: 

M5<\  and  kMS<\  (5.2) 

where  6  is  the  thickness  or  oscillatory  amplitude  of  the  airfoil,  whichever  is  the  larger;  and  k  is 
the  reduced  frequency  defined  as  k  =  . 
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According  to  his  large-Mach-number  expansion  theory,  Landahl  (Ref  5.13)  pointed  out  that 
Piston  Theory  amounts  to  ignoring  a  second  order  term  in  his  linear  amplitude  sequence.  Hence, 
the  valid  range  of  Mach  number  for  Piston  Theory  can  be  defined  by  the  criterion; 

S~^'^  <M <6~^  and  A/^»l  (5.3) 

In  terms  of  Tsien’s  Hypersonic  Similarity  parameter  (Ref  5.14)  AT,  where  K  =  MS ,  Eq  5.3 
reads: 


S^'^<K<\  (5.4) 

For  a  wedge  of  semi-angle  <7=10®,  K  falls  in  the  range  of  0.31  <  AT  <  1.0.  Inspection  of  results 
obtained  in  Fig  5.9  shows  that  the  valid  lower  bound  of  the  Mach  numbers  is  really  more 
restrictive  than  the  above  criterion  so  indicated,  whereas  the  upper  boimd  AT  <  1  is  less  restrictive. 


Fig  5.9  Surface  Pressure  of  a  Wedge  According  to  Various  Supersonic/Hypersonic  Models: 

T  =  tan  10®;  y  =  1.4 

While  the  condition  AT  <  1  for  Piston  Theory  may  be  somewhat  relaxed  to  include  regions  near  K 
«  1.0,  the  condition  kK  <  lof  Eq  5.2  is  a  rather  stringent  one.  For  unsteady  hypersonic  flow,  if 
K  «  0(1)  then  the  reduced  frequency  k  must  be  kept  very  small.  The  failure  of  Piston  Theory  in 
the  moderate  to  high  range  of  k  is  evidenced  by  the  Panel  Flutter  results  presented  in  the  work  of 
Chavez  and  Liu  (Ref  5.2). 

Following  the  suggestion  of  Morgen,  Herchel  and  Runyan  (Ref  5.15),  Rodden  and  Farkes  (Ref 
5.16)  have  arrived  at  a  generalized  expression  for  the  pressure  coefficients,  i.e.: 
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where  w  represents  the  piston  upwash. 
For  Piston  Theory: 


c,=l. 
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For  Van  Dyke’s  second  order  theory  (Ref  5.17): 


c,= 


m 


C2  — 


4^“ 


(5.6) 


(5.7) 


where  -  l,and  y  is  the  ratio  of  specific  heats. 

A  modified  Piston  Theory  is  recommended  to  replace  ci  and  C2  of  Eq  5.6  by  that  of  Eq  5.7 
rendering  an  extension  to  the  lower  Mach  number  region. 

The  Cl  and  cz  of  Van  Dyke  in  fact  were  first  obtained  by  Busematm  (Ref  5.18),  in  which  he  also 
included  a  third-order  term  based  on  a  consistent  expansion  of  the  simple  wave  theory,  i.e.: 

C3  =  +  c„  }  (5.8) 

6m 

where: 

=  y  +  1,  =  2r^  -ly  -  5,  c„  =  10{y  +  1) 

rf„=-12,  e„=8 

Following  Busemann,  Donov  (Refs  5.18,  5.19)  further  developed  a  comprehensive  theory  in 
which  he  obtains  series  expansion  solution  up  to  the  fourth-order  term  accounting  separately  for 
the  isentropic  part  and  the  rotational  part  due  to  simple  wave  and  shock  wave  respectively.  Here, 
Donov’s  third-order  term  including  shock  wave,  also  derived  independently  by  Carafbli  (Ref 
5.20),  reads: 


C3  — 


— +  cM"  +  dM^  +  e  } 
6Mm^  ^  ^ 


(5.9) 


where: 

-  Jr  +  if  .  3y=  -  12y  -  7  _  _  9(y  +  1) 

"  =  J.  * - i ^ - — 

d  =  -6,  e  =  4 
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In  passing,  it  is  noted  that  through  a  different  approach  Kahane  and  Lees  (Ref  5.21)  have 
obtained  a  correction  term  to  cs  of  Eq  5,8  resulting  in  essentially  the  same  C3  as  that  Eq  5.9. 
Therefore,  a  consistent  choice  of  Cp  would  be  to  adopt  Donov’s  series  and  Busemann’s  series  for 
flow  compression  and  expansion  respectively. 

In  the  analysis  that  follows,  we  remain  to  adopt  Lighthill’s  Piston  Theory,  Eq  5.6,  in  order  to 
simplify  the  present  approach. 

For  unsteady  flow  applications,  Eq  5.6  is  recast  into  the  form  of  pressure  differential  of  the  upper 
and  the  lower  wing  surfaces,  i.e.  ACp  =  Cp^  -  Cp^.  and  the  piston  velocity  w/C/„  is 

represented  by  two  terms,  i.e.  w  /  U„  =  where  denotes  the  thickness  distribution 

of  the  wing  and  w,  the  downwash.  Thus,  the  total  pressure  differential  ACp  can  be  expressed 
as; 

ACp  =  ACp_^  +  ACp  (5- 10) 

and  up  to  the  third-order  term: 

o  3 

-  -;4  £ 

and: 

ACp  =  %  f;  «c.  M"  (Aw.r’  +  Aw,'"  +4w,%  M  (5.10b) 


where: 
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For  non  lifting  airfoil  sections,  where  ),  =  (w^  )„. .  Eq  5. 10b  reduces  to  the  expression: 

ACp  =  -^[  (c,  )wi  +  ] 

M 


(5.11) 


Substituting  Eq  5.6  into  Eq  5.11  and  dropping  the  higher  order  terms  in  wi  yields  the  linear 
amplitude  version  of  Piston  Theory  as  a  special  case,  i.e.: 


ACp  = 


M 


2(/  +  l)w^  +  O'  +  OMw/ jwi 


(5.12) 
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We  shall  use  the  above  expression  to  develop  the  Unified  Lifting  Surface  Method. 
5.3  Hypersonic  Similarity  for  Thickness  Effects 
A  classical  Hypersonic  Similarity  (Ref  5.22)  can  be  expressed  as: 


where: 


/« 


r 


£_-l 

2 


is  the  universal  function  due  to  the  Prandtl-Meyer  expansion,  and: 


fn  =  K'^ 


fr+iV  1  p 

- -  +  — T 

I  4  j 


+ 


y +  1 

T~ 


1 


(5.13) 


(5.13a) 


(5.13b) 


is  the  universal  function  due  to  oblique  shock  waves,  where  K  =  MS  or  Mr . 

Clearly,  Eq  5.13a  is  the  basis  of  Lighthill’s  Piston  Theory  and  hence  of  Eq  5.5.  Eq  5.13b  was 
established  by  Tsien  (Ref  5.14)  and  Linnell  (Ref  5.23).  When  Eq  5.13b  is  expanded  up  to  the 
third-order  term,  the  coefficient  C3  corresponding  to  Eq  5.6  reads: 

c,  =  (5.14) 

^  32 


This  is  to  say  that  the  departure  between  Eqs  5. 13a  and  5. 13b  starts  firom  the  third-order  term  and 

(sy^  •  2y  "  5^ 

the  difference  of  which  amounts  to  Acj  =  ^  ^  <  0,  representing  the  difference  in 
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rotationality  due  to  shock  wave. 

It  is  desirable  to  extend  the  previous  third-order  theories  into  the  hypersonic  flow  regime  where 
K  >  0(1) .  Close  examination  of  them  reveals  that  the  Cp’s  of  these  third-order  theories  diverge 
drastically  as  K  increases  toward  the  Newtonian  limit  (Fig  5.9). 


Second-order  theories,  on  the  other  hand,  usually  result  in  one  half  the  value  of  Newtonian 
pressure,  whereas  Cp  of  Linear  theory  vanishes  at  the  Newtonian  limit. 

It  is  clear  that  Piston  Theory  does  not  yield  the  correct  limit  in  the  low  supersonic  end,  nor  does 
it  approach  the  Newtonian  limit  in  the  hypersonic  end.  Fig  5.9  shows  that  Piston  Theory  has  a 
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limited  valid  range  of  /n  AT  between  roughly  say  -1  to  at  most  0.5  (0.368  <K<  1.05),  for  a  wedge 
of  semi-angle  equal  to  10®. 

Clearly,  the  status  of  the  previous  third-order  theories  warrants  further  establishment  of  one  that 
is  imiformly  valid  and  covers  both  the  supersonic  and  hypersonic  limits.  In  the  present 
development  we  have  established  such  a  uniformly  valid  solution  by  means  of  a  strained 
parameter  technique  in  the  unified  supersonic-hypersonic  domain.  The  resulting  Cp’s  are  two 
composite  functions,  one  for  the  compression  waves  and  the  other  for  the  expansion  waves, 
which  can  be  generally  recast  into  a  pseudosimilar  form  as: 


(5.15) 


provided  that  the  coefficients  Ci,  C2  and  C3  could  be  suitably  chosen  fi-om  the  appropriate  third- 
order  theories. 


5.4  Unified  Supersonic/Hypersonic  Lifting  Surface  Method  of  ZONA7U 

The  matrix  equation  for  solving  ACp  on  wing-like  components  of  ZONA7U  can  be  expressed 
as: 


w,  =Dj.ACp^  (5.16) 

where  w,-  is  the  downwash  wingbox  due  to  structural  oscillation  {w,-}  =  F„ ,  F„  is  expressed 
in  Eq  3.26.  D,y  is  the  normal  velocity  influence  coefficients.  [D,y]=[MC]ww»  [MC]ww  is 
expressed  in  Eq  3.51  and  ACp  is  the  imsteady  pressme  difference  between  the  lower  and  upper 
surface  of  the /*  wingbox. 

It  has  been  shown  that  in  Ref  5.1  that  a  unified  NIC  matrix  can  be  constructed  by  superposing  a 
“nonlinear”  matrix  E,y  onto  the  matrix  Di,  based  on  the  principle  of  amplitude  perturbation,  i.e. 

w,  =[D,  (5.17) 

The  matrix  E,y  is  “nonlinear”  in  the  sense  that  it  contains  effects  due  to  nonlinear  functions  in 
wing  thickness  r  or  flow  incidence  These  nonlinear  effects  include  shock-induced 
rotationality  or  local  flow  expansions. 

There  are  a  number  of  approaches  that  could  provide  the  E,y  matrix.  For  example,  a  stripwise 
solution  could  be  provided  by  the  Perturbed  Euler  Characteristic  method  (Euler-Pec  method.  Ref 
5.2).  However,  the  simplest  and  most  expedient  approach  is  to  adopt  the  concept  of  Hayes- 
Lighthill’s  Piston  theory.  In  so  doing,  it  is  required  that  Ey  be  a  diagonal  matrix  whose  elements 
is  related  to  two  nonlinear  functions  in  Wo; 
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E/‘  =  f(w,-,r,M,K) 


(5.18) 


where  Wo  is  the  local  thickness  function  at  the  wingbox,  ^'is  the  ratio  of  specific  heat  of  gas 
and  K  (=Mt)  is  Tsien’s  hypersonic  similarity  parameter  (Ref  5.14).  The  functions /and  g  are 
uniformly-valid  solutions  in  the  unified  supersonic-hypersonic  domain  for  compression  and 
expansion  surfaces,  respectively.  By  means  of  strained  coordinate  technique  (Ref  5.24),  the 
function /is  derived  by  matching  the  Newtonian  impact  solution  with  then  Donov-Linell  series 
for  flow  compression  (Ref  5.19),  and  the  function  g  with  the  Busemann-Lighthill  series  for  flow 
expansion  (Ref  5.18).  Hence,  the  parameter  fiy  (Eq  5.17)  is  a  local  switching  operator  on  E,y  for 
suitable  adaptation  of  function  /or  g  on  each  panel.  It  should  be  noted  that  Eq  5.17  contains 
piston  theory  as  a  special  case,  in  which  D,y  =  0,  =  1  and  E,y  reduces  to: 

=  gK)  =  +  2(y+l)>v„  +  M(r  +  1)H'„^  (5.19) 

M 

As  commented  in  Ref  5. 1,  Piston  theory  inherits  two  undesirable  features.  First,  it  is  strictly  one¬ 
dimensional  model  which  provides  no  upstream  influence  whatsoever.  Second,  its  applicability 
in  Mach  number  range  is  ambiguous  for  it  does  not  approach  the  Ackeret  limit  in  the  low 
supersonic  end,  nor  does  it  approach  the  Newtonian  limit  in  the  hypersonic  end.  By  contrast,  the 
present  unified  solution  contains  both  limits.  In  the  Newtonian  limit,  where  M  approaches 
infinity  and  y  approaches  unity  simultaneously,  E,y  reduces,  as  expected,  to: 

E«‘‘  =  /(^o)  =  2^0^  and  20) 

E«"'  =  g(>^o)  =  0 

The  inadequacy  of  Piston  Theory  can  be  seen  in  the  case  of  an  oscillating  flap. 


Shock  Wave 


Fig  5.10  Oscillating  Leading-Edge  Flap  os  a  Thin  Wedge  Airfoil:  o  —  2° 

Fig  5.10  shows  an  oscillating  leading  edge  flap  of  thin  wedge  profile  (<r  =  2)  with  a  hinge  line 
located  at  the  quarter  chord.  Fig  5.11  shows  the  magnitude  and  phase  angle  of  the  unsteady 
pressure  along  the  profile  at  M=  5.0  and  k  =  0.5.  The  pressure  magnitude  on  the  flap  predicted 
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by  ZONA  51U  is  in  good  agreement  with  Euler-Pec  solution.  However,  the  phase  angle  of 
ZONA51U  follows  essentially  that  of  the  linear  theory,  which  disagrees  with  the  Euler-Pec 
solution.  This  is  expected  in  that  the  nonlinear  functions / and  g  of  Eq  5.18  only  contribute  to 
self-influenced  elements  of  E,y,  an  inherent  feature  of  Piston  theory  where  =  0,  it  essentially 
provides  no  phase  change;  hence,  its  prediction  is  inadequate  for  the  present  case. 
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Fig  5.11  Unsteady  Pressure  Distribution  for  an  Oscillating  Leading  Edge  Flap  with  Hinge  Line  at 

Quarter  Chord:  (Af  =  2.4,  k  =  0.5,  a  -  3®) 

5.5  AIC  Matrix  of  ZONA7U  for  Hypersonic  Wing-Body  Configuration 

The  inclusion  of  the  nonlinear  matrix  E;y  in  the  [MC]ww  matrix  suggests  that  the  total  normal 
velocity  influence  coefficient  matrix  for  wing-body  configuration  in  Eq  3.49a  can  be  modified 
as: 


[MCJbb  [MC]wb 

[mcJbw  [a^/C]^ +[//,•£,.] 

Since  solving  a  and  ACp  in  Eq  5.21  requires  the  matrix  inversion  of  the  modified  NIC,  the 

nonlinear  matrix  'Ey  indeed  gives  a  nonlinear  relationship  of  ACp  and  the  thickness  effect.  Also, 

the  procedures  fov  AIC  generation  described  in  Sec  3  can  be  directly  adopted  for  ZONA7U.  This 
leads  to  the  AIC  matrix  of  ZONA7U  to  be  in  the  same  form  of  that  of  ZONA7. 


AC. 


(5.21) 


95 


6.0  SPLINE  METHODS  FOR  SPLINE  MATRIX  GENERATION 


Aeroelastic  analysis,  as  an  interdisciplinary  problem,  requires  the  coupling  of  the  aerodynamic 
and  structural  responses.  In  practice,  the  requirements  to  generate  Ae  discretized  models  of 
these  disciplines  are  subject  to  different  engineering  considerations.  The  grid  of  the  discretized 
aerodynamic  model  is  usually  placed  on  the  external  surface,  whereas  that  of  the  structural 
model  is  placed  on  the  internal  load-carry  component.  This  gives  rise  to  the  data-transferred 
problem  between  two  computational  grid  systems.  This  would  amount  to  the  proper  transferal  of 
the  displacements  computed  in  the  structural  grid  to  the  aerodynamic  grid  and  that  of  loads  from 
the  aerodynamic  grid  to  the  structural  grid.  The  development  of  a  suitable  methodology  for 
solving  tWs  type  of  data  transferal  problem  is  by  no  means  a  trivial  task.  In  fact,  such  a 
methodology  should  be  further  developed  as  the  aerodynamic  and/or  structural  methods  advance. 

ZAERO  resolves  this  data  transferal  problem  by  mean  of  providing  a  spline  matrix  that  relates  or 
interpolates  the  displacements  at  the  structural  finite  element  grid  points  to  the  control  points  of 
aerodynamic  boxes.  This  spline  matrix  is  generated  by  a  spline  module  in  ZAERO  that  contains 
four  spline  methods  namely  the  beam  spline  method,  infinite  plate  spline  (IPS)  method,  thin- 
plate  spline  (TPS)  method  and  the  rigid-body  attachment  method.  These  four  methods  jointly 
assemble  the  total  spline  matrix  G  expressed  in  Eq  2. 15,  repeated  below: 

h  =  Gx 

where  h  is  the  “interpolated”  displacement  vector  at  aerodynamic  boxes,  including  the 
translational  displacements  and  their  slopes  with  respect  to  x.  Specifically,  h  is  in  the  order  of 
K-set  defined  in  Eq  3.60,  repeated  below: 


\\diere: 

i  represents  the  index  of  the  aerodynamic  box 
G  is  the  total  spline  matrix  relating  h  to  x 

X  is  the  displacement  vector  defined  at  the  structural  finite  element  grid  points 
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L  is  the  reference  length. 


Usually,  each  structural  finite  element  grid  point  has  six  degrees  of  freedom  (d.o.f.);  namely  T\, 
Ti  and  for  translational  displacement  as  well  as  i?i,  Ri  and  R^  for  rotational  displacement 
along  X,  y  and  2  directions.  Thus,  for  n  grid  points  there  are  6  x  n  degrees  of  freedom.  These  6  x 
n  degrees  of  freedom  are  defined  as  G-set  d.o.f.  that  can  be  expressed  as: 

li' 

T 
•‘2 

{x)G-s*t  =  n 

^2 

Us 

where  i  represents  the  index  of  the  structural  finite  element  grid  point. 

Once  the  spline  matrix  G  is  generated,  the  force  transferal  from  the  aerodynamic  boxes  to  the 
structural  finite  element  grid  points  (from  K-set  to  G-set)  can  be  achieved  by  the  transpose  of 
matrix  G.  This  has  been  shown  in  Eq  2.16,  repeated  below: 

{FJo-./  =  [Gf  {V^K-se. 

In  this  section,  we  will  discuss  the  theoretical  derivation  of: 

•  Infinite  Plate  Spline  (IPS)  method 

•  Thin-Plate  Spline  (IPS)  method 

•  Beam  Spline  method 

•  Rigid-Body  Attachment  (RBA)  method 

•  Matrix  Assembly  of  the  Total  Spline  Matrix 

6.1  The  Infinite  Plate  Spline  (IPS)  Method 

The  IPS  method  was  first  proposed  by  Harder  and  Desmarais  (Ref  6.1),  which  was  a  significant 
improvement  over  the  two-dimensional  (2D)  interpolation  method  of  Rodden  (Ref  6.2).  This 
development  was  motivated  by  the  advent  of  lifting  surface  methods  in  aerodynamics  at  that 
time,  which  required  a  2D  interpolation  method  such  as  IPS.  The  2D  surface  is  defined  as  the 
plane  of  the  lifting  surface.  Therefore,  IPS  is  ideally  suited  for  displacements  and  forces 
transferal  of  wing-like  components.  Today,  EPS  is  one  of  the  more  popular  method  of 
interpolation  used  in  aerospace  industry. 
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Consider  a  set  of  discrete  structural  grid  points  (;c/,  for  i=  \,2,  ...,N  lying  within  a  2D 
domain  with  Cartesian  coordinates  x  and  y.  Each  grid  point  defines  the  vertical  position 
coordinate  of  the  “deformed”  surface.  IPS  solves  the  partial  differential  equation  of  equilibrium 
for  an  infinite  plate  with  uniform  thickness.  The  deformation  of  the  infinite  plate  satisfies  the 
given  deflection  w,{xi,  yt)  at  the  N  structural  grid  points.  Once  the  partial  differential  equation  is 
solved,  the  deflection  at  other  points,  for  instance  the  aerodynamic  points,  on  the  plate  can  be 
determined. 


The  governing  equation  of  an  infinite  plate  with  bending  stiffness  reads: 
DV^W  =  q 


(6.1) 


where  W  is  the  plate  deflection,  D  is  the  plate  bending  rigidity,  and  q  is  the  distributed  load  on 
the  plate.  Introducing  polar  coordinates,  x  =  rcosO ,  y  =  r sin ^ ,  so  that  in  polar 
coordinates  is  given  by: 


r  dr 


(6.2) 


and  considering  the  deflection  due  to  a  point  load  P  at  the  origin  of  the  coordinate  system,  a 
solution  of  Eq  6. 1  can  be  written  as: 


W(r)  =  A  +  Br^ 


l67rDj 


Inr^ 


(6.3) 


where  A  and  B  are  undetermined  coefficients. 

For  N  point  loads  at  the  given  location  (x,,  y,),  for  /  =  1,  2,  ...,  AT  in  the  2D  space,  the  total 
deflection  can  be  obtained  by  superimposing  the  fundamental  solution  of  Eq  7.1  such  that: 


W{x,y)  =  Y,  [a  +  a  +  Pi  ^i  Inr/  ) 


(6.4) 


where: 

p 

Ai,  Bi  and  F,  =  — ■ —  are  undetermined  coefficients, 

'  16^D 

and: 

+  (y-y,)^ 

For  the  purpose  of  determining  these  undetermined  coefficients  one  needs  to  use  certain 
information  about  the  solution.  Harder  and  Desmarais  showed  that  by  expanding  Eq  6.4  for 
large  values  of  r,  one  obtains  terms  of  order  r\  r,  1,  Mr,  etc.,  along  with  terms  of  order  In  r 
In  y.  In  y,  etc.: 
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(6.5) 


fV(r,0)  =  rMnr^  YF,  +  YB,  -  2rlnr^  cosB+;;,  sinB)Fi 
/=!  <=1  /=! 
iV  , 

-  2r  y(.,  cosB  +  ^i  sinB)(Fj +B,)  +  +3;,  )F,  +  ... 

M  <=1 

For  removing  the  singularity  at  r  =  00 ,  coefficients  of  the  terms  of  order  In  and 

r  In  7^  must  vanish.  This  gives; 

=  0 

fal 


(6.6) 


N 


(6.7) 

'Ey,  F,  =  0 

1=1 

(6.8) 

N 

Z«,  =  0 

1=1 

(6.9) 

Here  Eq  6.6  can  be  recognized  as  the  discrete  force  equilibrium  equation,  which  eliminates  terms 
of  order  In  r^.  Eqs  6.7  and  6.8  are  discrete  moment  equilibrium  equations  and  eliminate  terms 
of  order  r  In  Finally,  Eq  6.9,  the  physical  significance  of  which  is  not  clear,  serves  to 
eliminate  terms  of  order  r^. 

Eqs  6.6  through  6.9  result  in  linear  deflection  at  infinity.  For  extrapolation,  this  implies  that 
linear  deflection  of  the  aerodynamic  points  occurs  only  if  they  are  located  far  from  the  domain  of 
the  structural  grid  points.  A  solution  to  the  general  spline  problem,  formed  by  superimposing 
solutions  of  Eq  6.1,  is  given  by: 


Hx,y)  =  ^0 


N 

+  X  +  y  +  2i:,(x,y)F, 


1=1 


(6.10) 


where: 

K,{x,y)  =  r/  Inr,^ 

=  (x-jc,)2  + 

where  flo,  oi  and  02  are  xmknowns  given  by: 
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(6.11) 

i=i 

N 

/»! 

(6.12) 

N 

a,  =-2£B,y, 

i=l 

(6.13) 

Note  that  the  //+  3  unknowns  in  Eq  6.10  can  be  determined  from  application  of  side  conditions 
found  in  Eqs  6.6  -  6.8  along  with  setting  the  deflection  at  the  point  to  its  known  value  W}. 
Viz., 

W,  =  +  a^  X,  +  02  yt  +  Fj  for  i  =  1, 2, V  (6.14) 


where: 


K.. 


(6.15) 


r,/  =  +  (y,-yjf 

is  the  square  of  the  distance  between  known  points  (xi,  yi)  and  (xj,yJ). 

Eq  6.14  and  the  side  conditions  foimd  in  Eq  6.6  -  6.8  can  now  be  expressed  in  matrix  form  as: 


{W}  =  [R]{a}  +  [K,]{F} 


(6.16) 


and: 

[Rf  {F}  =  0 


where: 


'  F,' 

{W)  =  ^ 

W2 

(F}  =  - 

V _ 

r 

’^1' 

yi' 

{a}  =  < 

[R]  = 

< 

1 

►  ^ 

X2 

►  * 

y2 

► 

1 

k  j 

/n. 

.>'n. 

(6.17) 


(6.18) 
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Thus,  the  deflection  JV/,  for  i  =  1, 2,  ...,N  can  be  detennined  using  Eqs  6.16  and  6.17  along  with 
the  above  definitions.  Combining  Eqs  6. 16  and  6. 17  gives: 


/ - 

o 

I - 

t-H 

o 

o 

o 

1 _ 

r 

0 

0  0  0  .Xj  •••  Xj, 

ai 

0 

0  0  0  y,  —  y^ 

a2 

w^ 

.  = 

1  x^  yi  0 

W2 

1  ^2  yi  ^21  *'■  ^2N 

••• 

A 

_  i  .Vn  ^N1  *"  ®  J 

=  [C][P] 


(6.19) 


The  interpolation  to  any  point  in  the  2D  plane  is  then  achieved  by  evaluating  w(x,  y)  from  Eq 
6.10  at  the  desired  points.  Thus,  for  a  given  aerodynamic  point  its  displacement  is: 


O' 

0 

0 


h(jCt,yjt)  (1,  •^ic,  y^f 


k,2> 


)  [C] 


-1 


w, 


and  its  streamwise  slopes  is: 


O' 

0 

0 


=  (0,  -1,  0,  DATj^i,  dk:^2>  •••>  )  [CJ 


-1 


w, 


w. 


(6.20) 


(6.21) 


where: 

DK^,  =  =  -2(^1  -  +  ’"'■k,.’) 

aij 

Eqs  6.20  and  6.21  can  be  rewritten  as: 
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(6.22) 


=  [Skl  i 


w, 


w, 


yv, 


nj 


and: 


hXx,,y,)  =  [DsJ 


w. 


w. 


(6.23) 


where  [s^]  and  [Z)s,j]  are  resulted  from  the  matrix  multiplications  in  Eqs  6.22  and  6.23, 
respectively,  but  with  the  first  three  columns  of  [C]"'  being  removed. 

The  solutions  of  h(x,j,j'k)  and  h'(.:Cic,)'k)  exist  only  if  the  matrix  C  is  non-singular.  The 
singularity  in  matrix  C  occurs  when: 

•  all  structural  grid  points  are  aligned  along  a  line.  This  is  obvious  since  a  line  fails  to  define  a 
plane. 

•  two  or  more  than  two  structural  grid  points  have  the  same  x  and  y  location. 

To  perform  the  IPS  method,  it  is  required  that  all  structural  grid  points  and  aerodynamic  points 
are  located  on  the  same  plane.  This  plane  is  called  “spline  plane”.  Normally,  the  plane  of  the 
lifting  surface  (or  the  mean  plane  of  the  wing-like  component)  is  selected  as  the  spline  plane. 
However,  for  structural  grid  points  located  in  3D  space,  these  structural  grid  points  may  not  be 
necessarily  located  on  the  plane.  In  this  case,  it  is  required  to  project  the  structural  grid  points 
onto  the  spline  plane  along  the  normal  direction  of  the  spline  plane.  This  can  be  done  by 
transforming  the  structural  grid  point  locations  to  a  local  coordinates  whose  x-y  plane  coincide 
with  the  spline  plane. 

Singularity  in  matrix  C  appears  in  case  two  structural  grid  points  shear  the  same  x  and  y 
locations  on  their  projected  position,  even  if  their  original  positions  in  the  3D  space  are  not  the 
same.  Therefore,  in  this  kind  of  situation,  one  of  the  two  grid  points  must  be  excluded  from  the 
selection  of  the  structural  grid  points. 

Another  important  concept  needed  to  be  addressed  here  is  that  the  IPS  method  is  a  scalar 
operator.  This  is  to  say  that  for  a  given  set  of  normal  displacements  at  structural  grid  points  the 
IPS  method  results  the  displacements  at  the  aerodynamic  points  only  along  the  normal  direction 
of  the  spline  plane.  For  instance,  the  deflection  W  in  Eq  6. 1  represents  the  normal  displacement 
at  structural  grid  point.  Therefore,  the  deflection  h  in  Eq  6.20  at  aerodynamic  point  is  also  the 
normal  displacement.  However,  for  the  deflection  along  other  directions  one  finds  out  that  the 
IPS  method  can  be  applied  in  the  same  way  as  that  of  the  normal  displacements.  Let  u  be  the 
streamwise  deflections  and  v  be  the  lateral  deflections,  Eq  6. 1  can  be  rewritten  as: 
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DV*u  =  q 
DVS  =  q 


(6.24) 


(6.25) 


This  time,  D  is  not  the  plate  bending  rigidity  but  represents  some  in-plane  flexural  rigidity. 
Solving  Eqs  6.24  and  6.25  results  the  solutions  that  are  identical  to  Eqs  6.22  and  6.23  but  with  W 
replaced  by  u  for  streamwise  displacement  and  by  v  for  lateral  displacement.  This  also  indicates 
that  for  6  d.o.f  at  each  structural  finite  element  grid  points,  namely,  Ti,  T2,  T3  of  translational 
displacements  and  Ru  Rz,  R3  for  rotational  displacements,  Eqs  6.22  and  6.23  can  be  generalized 
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(6.26) 


The  zero  matrices  associated  with  (Ri) ,  {Rj}  and  {R3}  in  the  right  hand  side  of  Eq  6.26  are 
resulted  from  the  IPS  formulation  which  does  not  involve  the  structural  rotational  displacements. 


Applying  Eq  6.26  for  M aerodynamic  points  by  letting  the  index  k  be  k  =  1, 2, . . . ,  Myields: 


{h}  =  [G^]{x} 


where; 


k=l,2. 


M 


(6.27) 


(6.28) 
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(6.29) 


(6.30) 


Thus,  [G  IPS  ]  represents  an  EPS  method  generated  spline  matrix  that  interpolates  the 

displacements  of  the  6  d.o.f  at  N  structural  grid  points  to  the  displacements  along  x,  y,  z 
directions  as  well  as  their  streamwise  slopes  at  A/ aerodynamic  points. 

Finally,  some  important  remarks  about  the  EPS  method  are  summarized  below; 


•  The  spline  matrix  becomes  singular  when  two  or  more  than  two  structural  grid  points  shear 
the  same  x  and  3^  location  of  their  projection  position  on  the  spline  plane. 

•  The  spline  matrix  is  singular  if  all  structural  grid  points  are  aligned  along  a  line. 

•  Linear  extrapolation  occurs  only  if  the  aerodynamic  points  are  far  from  the  domain  of  the 
structural  points.  Otherwise,  distortions  or  oscillations  may  appear  in  the  extrapolated 
regions. 

•  The  PS  method  is  a  scalar  operator.  For  a  given  set  of  displacements  along  one  direction, 
the  PS  method  does  not  recover  displacements  along  other  directions. 

6.2  The  Thin-Plate  Spline  (TPS)  Method 

The  TPS  method  (Ref  6.3)  is  a  generalization  of  the  PS  method  by  incorporating  some  three- 
dimensional  aspects.  TPS  provides  a  means  to  characterize  an  irregular  surface  by  using 
functions  that  minimize  an  energy  functional.  The  derivation  is  entirely  analogous  with  the  PS 
method  with  the  addition  of  the  third  coordinate.  Eq  6.4  becomes: 

W(_x,y,z)  =  f  {a,+B,  +  F,  Inr/  )  (6.31) 

i=l 


where: 
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r/  =  +  iy-y,f  + 


The  boimdary  conditions  at  infinity  are  similar  to  these  of  IPS  but  with  the  addition  moments  in 
the  third  axis: 

Sa  =  =  Z^.  =  'Zy<  F,  =  -p;  («-32) 

Expanding  Eq  6.31  for  large  r,  terms  of  order  rlnr,r  and  r  In  r'  can  be  eliminated  by  applying 
the  boundary  conditions  expressed  in  Eq  6.32.  This  leads  to: 

JV 

w(x,y,z)  =  +  a^x  +  a2y  +  a^z  +Y^^i(.^>y>^)Fi  (6.33) 

i=i 


where: 

K,(x,y)  =  Inr,^ 

In  order  to  solve  the  unknown  coefficients,  one  can  introduce  a  matrix  notation  such  that  Eqs 
6.32  and  6.33  become: 


{W}  =  [R]{a}  +  [K,]{F} 
[Rf  (F)  =  0 

where: 

K^.  =  r/  Inr/ 


and: 

r/  =  (xj-x^f  +  (yj-yif  +  {zj-Zif 
The  matrix  R  is: 


[R]  = 


1 

<•  ^ 
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►  < 

►  ■< 

72 

►  ^ 

^2 

► 

1 

.7n. 

{W}  = 


'W,  ' 

^7 

«2 

.^3. 

F7 

2 

< 

(a)  =  ^ 

(F)  =  . 

2 

(6.34) 


(6.35) 


(6.36) 
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For  a  given  coordinates  of  the  N  structural  grid  points,  Eq  3.34  becomes: 


[0]  [R]" 
[R]  [K,]j 


{IS}  ■  '"{SI 


(6.37) 


The  interpolation  to  any  point  in  the  3D  space  is  then  achieved  by  evaluating  w^x,  y,  z)  of  Eq 
6.33  at  the  desired  point.  Thus,  for  a  given  aerodynamic  point  {x^ ,  A ,  ) ,  its  displacement  is: 


~  (^>  -^k"  A*  ^k»  •^k,l>  ^k,2>  -^kn  )  [^J 


\  (0} 

.W. 


and  its  streamwise  slopes  is: 


hK>A,^k)  =  (0.  -1.  0,  0,  [C] 


•{ 


(0)1 

{W)J 


where: 


DK^,  =  =  -2{x,  -  .x,)(l  +  lnr,,,=^) 


Eqs  6.38  and  6.39  can  be  rewritten  as: 

=  (Al  {W} 

and: 

h'(A:^,A.z,)  =  [DsJ  {W} 


(6.38) 

(6.39) 

(6.40) 

(6.41) 


where  [s^]  and  [Ds,j]  are  resulted  from  the  matrix  multiplication  in  Eqs  6.38  and  6.39, 
respectively,  but  with  the  first  four  columns  of  [C]"‘  being  removed. 

Unlike  the  ff  S  method,  TPS  does  not  require  a  spline  plane.  Therefore,  TPS  can  be  considered 
as  a  3D  spline  method  that  performs  the  interpolation  based  on  a  set  of  structural  points  located 
in  a  3D  space. 

Similar  to  IPS,  the  TPS  method  has  several  restrictions: 

•  The  matrix  [C]  in  Eq  6.37  is  singular  if  two  structural  grid  points  have  the  same  x,  y  and  z 
locations. 

•  The  matrix  [C]  is  singular  if  all  structural  grid  points  are  located  on  the  same  plane.  In  this 
case,  the  sub-matrix  [R]  becomes  singular  that  leads  to  a  singular  matrix  of  [C]. 
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TPS  is  also  a  scalar  operator.  Similar  to  IPS,  Eq  6.41  can  be  generalized  into  a  spline  matrix 
such  that: 

{11}  =  [GtkK*} 

where  {h}  and  {x}are  the  same  as  Eqs  6.28  and  6.29,  respectively,  and  [Gtps]  has  the  same  form 
as  Eq  6.30. 

Since  there  is  no  requirement  of  spline  plane,  for  interpolation  of  displacements  at  3D  structural 
grid  points,  TPS  outperforms  the  PS  method.  However,  if  all  structural  grid  points  are  on  the 
same  plane,  the  TPS  method  breaks  down.  In  theory,  for  this  kind  of  situation  the  formulation  of 
TPS  can  be  reduced  to  that  of  PS.  But  such  a  reduction  is  not  straight  forward.  Therefore,  in 
practice,  the  PS  method  should  not  be  treated  as  a  special  case  of  TPS.  Thus,  for  all  structural 
grid  point  on  the  same  plane,  the  use  of  that  PS  method  is  suggested. 

6.3  The  Beam  Spline  Method 

It  is  often  that  a  high  aspect  ratio  wing  structure  is  modeled  by  a  beam-type  finite  elements  along 
the  elastic  axis  of  the  wing,  or  a  body  by  beam-type  finite  elements  along  the  center  line  of  the 
body.  On  the  other  hand,  the  aerodynamic  boxes  of  the  wing-like  and  body-like  components  in 
the  aerodynamic  model  are  generally  located  in  a  3D  space.  The  Beam  Spline  Method  is 
designed  to  specially  handle  this  kind  of  spline  problem. 

The  Beam  Spline  Method  solves  the  partial  differential  equation  of  equilibrium  for  an  infinite 
beam  with  uniform  bending  and  torsion  stiffness.  For  bending  deflections,  it  satisfies  the  given 
displacements  and  slopes  at  the  N  structural  grid  points.  For  torsion,  it  satisfies  the  given  twists. 
The  N  structural  grid  points  are  assumed  to  be  located  along  a  line  called  the  “spline  axis”.  In 
the  present  formulation  of  the  Beam  Spline  Method,  the  spline  axis  is  defined  as  the  y  axis  of  a 
user  specified  local  coordinate  system. 

After  the  unknowns  of  the  infinite  beam  equation  are  determined,  the  displacements  and  slopes 
at  the  aerodynamic  points  can  be  obtained  by  rigid  connections  between  the  beam  and  the 
aerodynamic  points  (see  Fig  6.1).  To  derive  the  solution  of  the  infinite  beam  equation,  it  is 
required  first  to  transform  all  structural  grid  points  and  aerodynamic  points  into  the  spline  axis 
coordinate  system.  This  can  be  done  by  computing  a  transformation  matrix  [Ts]  such  that: 


f  > 

X 

^y 

^  =  [TJ  - 

z' 

V.  J 

2 

(6.43) 


where: 

x' ,  y  and  z'  are  a  local  coordinate  system  whose  y'  axis  is  the  spline  line 


and: 
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X  ,  y  and  z  are  the  basic  coordinate  system  whose  x  axis  is  parallel  to  the  streamwise 
direction. 


z 


y*  (spline  line) 


z 


X 


(b) 


F%6.1  Spline  Axis  Coordinate  System  (a)  Spline  Axis  Along  the  Elastic  Axis  of  Wing-Like 
Component  (b)  Spline  Axis  Along  the  Center  Line  of  Body-Like  Component 
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Bendim  Stiffness  About  the  x  *  Axis 

The  partial  equation  of  an  infinite  beam  with  constant  bending  stiffness  about  x'  axis  reads: 


El 


d^W 

d/ 


(6.44) 


where: 

W  is  the  beam  deflection  parallel  to  the  z'  axis 
E I  represents  the  bending  stiffiiess 
g  is  a  distributed  transverse  load 

and: 

Mis  a  distributed  moment 

For  N  structural  grid  points  located  along  the  beam,  the  solution  is  found  by  superimposing  the 
fundamental  solutions: 
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(6.45) 

(6.46) 


Applying  Eqs  6.45  and  6.46  for  the  given  W  and  0^  at  N  structural  grid  points  and  imposing  the 
boimdary  condition  at  infinity  for  linear  function  ofW: 
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For  a  given  aerodynamic  point  y^,  Zy.)  exactly  located  at  the  beam  so  that  xi-  —  Zk  =  0,  the 
interpolated  displacement  W(0,>^k,0)  can  be  obtained  by  evaluating  W(y')  of  Eq  6.45  at  the 
desired  point: 
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W(0.y„0)  = 


<  ,  l^kll  kkkl  l?w|  |?kll 

’  -^k’  i-.r»  fr’  ’  n  irr  ’  /i  rr  ’ 


12  FI  12  El 


7k2  kiczl  ^ky  b  kv| 


4£'/ 


4  £7 


12£7’  4£:/ 

{0> 


'  [C]-'  j  {W} 

[{^k}J 


J^(0.n.0)  = 


J  Vu  l^kll  ^k2  |7k2|  Vw  l^kyj  l^kll 

’  4£7  ’  AEI  A  El  ’  2£7’ 


|^k2|  ^  _  bky 


2£:/’  ’  2EI 


'  {0}  ' 

[C]-‘  - 

{W}  ► 

ttCO^a.o)  = 


'o.o.  W 


k2 


sign(7k2) 

2  El 


2Er  2EV  ’  2EV 

sign(7kAr) 


W  sign(77fct) 

2E1  ' 


2  El 


'  [C]-‘ « 

■  {0}  ' 

{W}  ► 

J 

.  (^k) . 

Eqs  6.50, 6.51  and  6.52  can  be  simplified  and  rewritten  as: 

w(o.x.o)  =  k  »k]l' 

^(O.X.O)  =  [Z)s»  DRjj 
>’k  I 


{W} 


<m 


{W) 


|jM,0)  =  [csk 


where: 


(6.50) 


(6.51) 


(6.52) 


(6.53) 


(6.54) 


(6.55) 


[s,  ],  [Z)s,  DR^  ].  [Cs.  CRt  ]  are  resulted  from  the  matrix  multiplications 

in  Eqs  6.50, 6.51  and  6.52,  respectively,  but  with  the  first  two  columns  of  [C]'^  being  removed. 


Bendins  Stiffness  About  the  z  *  Axis 

The  partial  equation  of  an  infinite  beam  with  constant  bending  stiffness  about  the  z'  axis  reads: 
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(6.56) 


d*V  dM 


where  U  is  the  beam  deflection  parallel  to  the  x'  axis. 


The  equivalence  of  Eq  6.56  and  Eq  6.44  suggests  that  the  solution  of  Eq  6.44  can  be  adopted  for 
Eq  6.56.  Thus,  Eqs  6.53, 6.54  and  6.55  can  be  rewritten  as  the  solutions  for  U: 

{U}1 

ml 


A  MXXW  V/ .  ft 

u(o,y;.o)  =  [sj  Rjj 
r(O.X.O)  =  [Ds,  DRjj 


dy\ 


{U} 


0(O.X.O)  =  [Cs.  cr,]{2} 

Torsion  Stiffness  About  the  v  *  Axis 

The  equation  of  infinite  beam  with  constant  torsion  stiffiiess  reads: 


GJ 


=  -My 


where: 


0  is  the  twist  about  y’ 


My  is  a  distributed  torque. 


For  N  structural  grid  points  located  along  the  beam,  the  solution  of  Eq  6.60  is: 

w  = 


'  \y{-y\  \y2-y\  ^ 

_  ’  2GJ  ’  2GJ  2GJ 


(6.57) 


(6.58) 


(6.59) 


(6.60) 


(6.61) 


The  equilibrium  condition  requires: 

^My.  =  0  (6.62) 

Applying  Eq  6.61  for  the  twists  at  N  structural  grid  points  and  imposing  the  equilibrium 
condition  yield: 
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{M.} 


(6.63) 


(M.) 


For  a  given  aerodynamic  point  (0,  y'^ ,  0)  at  the  beam,  the  interpolated  twist  is: 


0yi0,yl,0)  =  [TJ  {0^} 


(6.64) 


where: 


[TJ  is  resulted  from  tire  combination  of  Eq  6.61  and  Eq  6.63. 


Displacement  of  Aerodynamic  Point  in  3D  Space 


The  displacement  of  a  given  aerodynamic  point  in  3D  space  is  obtained  by  assuming 

a  rigid  bar  connecting  the  point  and  the  beam.  This  leads  to: 


vl  =  [Gi(Sk,Rk,DSk,DRk,T,^,Xfc,zj]{x'} 


(6.65) 


=  [G,(£>s.,i>R,.Tjl{x'} 


(6.66) 


du 

dv 

dw 


=  [G3(Z)Sk,DRk,CSk,CRk,Xt,zj]  {x'} 


(6.67) 
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=  [G,(Osj,i)Rj.T„)]{»'} 


(6.68) 


du 


dv 
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where: 
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{w} 

R) 

{0,) 

.  {^z)  . 


(6.69) 


Gi,  G2,  G3  and  G4  can  be  derived  from  Eqs  6.53  to  6.55,  Eqs  6.57  to  6.59,  and  6.60  as 
well  as  with  a  rigid  bar  connecting  the  aerodynamic  point  (;c^ ,  )  to  the  beam. 


Transformation  from  Spline  Axis  Coordinate  to  Basic  Coordinates 

Eqs  6.65  to  6.69  are  derived  in  the  spline  axis  coordinates.  The  transformation  matrix  [Ts] 
relating  the  spline  axis  coordinates  to  die  basic  coordinates  is  defined  in  Eq  6.43.  Therefore,  Eq 
6.69  can  be  transformed  to  Ti,  T^,  Ts  and  R\,  R2,  R3  of  the  given  displacements  defined  in  the 
basic  coordinates  by: 


(x'} 


[TJ  0 
0  [TJ 


{X} 


(6.70) 


where  (x)  is  expressed  in  Eq  6.29. 

The  interpolated  displacement  u,  v  and  w  can  be  transformed  to  h^,  hy  and  hz  by: 


U 

•  =  [T.r  < 

V 

A. 

k 

w 

Substituting  Eqs  6.70  and  6.71  into  Eq  6.65  yields: 


(6.71) 
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(6.72) 


=  [Tsf  [G.] 
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The  streamwise  slopes  can  be  derived  from: 
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(6.73) 


Substituting  Eqs  6.66  to  6.69  and  combining  the  resultant  equation  with  Eq  6.72  yield  a  spline 
matrix  relating  the  6  d.o.f.  displacements  at  N  structural  grid  points.  Repeating  this  process  for 
A/ aerodynamic  points  by  letting  the  index  k  vary  from  1  to  M gives: 

=  (6.74) 


where: 

(h)  and  {x}  are  defined  in  Eqs  6.28  and  6.29,  respectively. 


and: 

[GBeaJ  is  a  6  X  M  spline  matrix  generated  by  the  beam  spline  method. 

Some  important  remarks  about  the  Beam  Spline  Method  are  discussed  as  follows: 

•  One  of  the  basic  assumptions  of  the  Beam  Spline  Method  is  that  all  structural  grid  points  are 
located  along  the  spline  axis.  Errors  may  be  introduced  if  some  of  the  grid  points  are  off  the 
axis. 

•  Similar  to  the  IPS  and  TPS  methods,  the  spline  matrix  becomes  singular  if  two  structural  grid 
points  are  coincided  at  their  projected  position  at  the  spline  axis. 

•  Linear  extrapolation  occurs  only  if  the  aerodynamic  points  are  far  away  from  the  domain  of 
the  structural  points.  Otherwise,  distortions  or  oscillations  may  appear  in  the  extrapolation 
regions. 

6.4  The  Rigid-Body  Attachment  (RBA)  Method 

The  Rigid-Body  Attachment  (RBA)  method  is  developed  for  the  cases  where  no  structural  model 
exists  for  a  particular  aerodynamic  component.  For  example,  an  underwing  store  may  be 
modeled  as  a  point  mass  in  the  structxiral  finite  element  model,  but  its  surface  may  be  represented 
by  a  detailed  aerodynamic  panel  model.  Since  no  information  of  the  structural  connectivity 
between  the  point  mass  and  the  store  surface  is  given,  it  is  assumed  that  all  aerodynamic  points 
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are  attached  to  the  point  mass  by  a  rigid  body.  For  a  given  aerodynamic  point  located  at  x^,  yk, 
Zk,  its  displacements  can  be  expressed  as: 
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and  its  streamwise  slopes  are: 
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(6.75) 


(6.76) 


where  x,  y  and  z  are  the  location  of  the  structural  grid  point  to  which  the  aerodynamic  points  are 
attached. 


Applying Eqs  6.75  and  6.76  for  M aerodynamic  points  by  letting  the  index  khtk-  1,2,  ...,M 
yields: 


{M  =  [Grba]{x} 


(6.77) 


where: 

{h}  and  {x}  are  defined  in  Eqs  6.28  and  6.29,  respectively. 


and: 

[Grba]  is  a  6  X  Mby  6  spline  matrix  generated  by  the  Rigid-Body  Attachment  method. 

Since  there  is  no  matrix  inversion  involved  in  the  Rigid-Body  Attachment  method,  [Grb^] 
cannot  be  singular.  Unsatisfactory  results  of  the  Rigid-Body  Attachment  method  is  usually 
caused  by  the  over-simplified  finite  element  model.  For  instance,  using  the  transposed  of 
[Grba]  force  transferal  from  aerodynamic  points  to  the  structural  point  implies  that  the 

entire  aerodynamic  forces  on  the  aerodynamic  component  are  lumped  at  a  single  structural  point. 
This  will  create  a  highly  concentrated  load  on  the  stmcture  which  may  not  be  realistic.  To  avoid 
this  kind  of  problem,  it  is  suggested  to  refine  the  structural  finite  element  model  by  introducing 
more  grid  points  and  using  other  spline  methods  such  as  TPS. 
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6.5  Matrix  Assembly  of  the  Total  Spline  Matrix 

The  generation  of  the  total  spline  matrix  that  relates  the  displacements  at  the  structural  G-set 
d.o.f  to  the  aerodynamic  K-set  d.o.f.  is  performed  on  a  component-by-component  basis.  The 
selection  of  the  spline  methods  depends  on  the  type  of  the  component  in  the  aerodynamic  model 
(i.e.  wing-like  or  body-like  components)  and  the  type  of  elements  (beam  or  plate  element)  in  the 
structural  finite  element  model.  Within  each  component,  the  use  of  several  splines  is  also 
allowed.  For  example,  a  model  may  use  the  beam  spline  method  for  the  inboard  section,  the  TPS 
method  for  the  outboard  section  and  the  IPS  method  for  the  aileron  of  a  wing-like  component. 
Separation  into  sub-regions  of  a  component  allows  discontinuous  slopes  (e.g.  along  the  wing- 
aileron  hingeline)  and  discontinuous  displacements  (e.g.  along  the  inboard  and  outboard  edges  of 
ailerons).  Therefore,  the  spline  matrices  generated  by  different  spline  methods  on  different  sub- 
regions  are  the  sub-matrices  of  the  total  spline  matrix.  Eq  6.78  shows  an  example  of  such  sub¬ 
matrix  arrangement: 

[G^s]  -  0  •••  0 

•  *  • 

«  •  * 

•  •  ♦ 

0  •••  [Cr-j-pg]  •••  0 

•  •  * 

•  •  ♦ 

•  •  ♦ 

[G,ps]  -  0  •••  0 

0  •••  0  •••  [Grba] 

The  degrees  of  freedom  in  (h)  and  {x}  are  grouped  based  on  the  aerodynamic  and  structural 
points,  respectively,  involved  in  each  spline.  Therefore,  it  is  required  to  rearrange  the  rows  and 
columns  of  the  total  spline  matrix  in  Eq  6.78  according  to  the  K-set  and  G-set  degrees  of 
fi-eedom  of  the  aerodynamic  model  and  structural  finite  element  model,  respectively. 

Finally,  it  should  be  noted  that  some  of  the  degrees  of  freedom  in  the  structural  finite  element 
model  could  be  specified  in  local  coordinate  systems.  It  is  these  coordinate  systems  that  define 
the  displacement  vectors  computed  by  the  structural  finite  element  method.  Therefore,  a  final 
transformation  of  (x)  for  these  degrees  of  freedom  is  required.  Let  [Tm]  be  a  transformation 
matrix  that  relates  the  local  to  the  basic  coordinate  systems,  then  the  spline  matrix  becomes  the 
total  such  that: 

=  IG]  Wo-,. 

where: 

tG]  =  [G]  [t„] 

X  is  the  structural  displacements  defined 

{>■)  =  [tJ  W 

and: 

[G]  is  the  total  spline  matrix  for  structural  displacement  defined  in  basic  coordinates. 
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7.0  FLUTTER  SOLUTION  TECHNIQUES 


ZAERO  provides  flutter  solutions  for  aeroelastic  analysis.  Several  flutter  methods  are  included 
in  ZAERO.  In  this  section,  we  will  discuss  the  solution  technique  of  various  flutter  methods. 
These  methods  include: 


•  the  K-method 

•  the  P-K  method 

•  the  g-method 

Technical  merits  and  theoretical  validity  of  each  method  will  also  be  discussed. 

In  Sec  2,  we  have  derived  the  flutter  matrix  equation  in  the  Laplace  domain  (Eq  2.9)  in  terms  of 
the  generalized  mass  matrix,  M,  generalized  stiffness  matrix  K,  and  the  generalized  aerodynamic 
force  matrix  Q.  Eq  2.9  reads,  repeated  below; 


q  =  0 


(7.1) 


Introducing  a  non-dimensional  Laplace  parameter p  defined  in  Eq  2.23,  repeated  below: 

P  =  ^  =  0*  +  fi:)  (7.2) 

where  k  is  the  reduced  frequency  defined  in  Eq  2. 13,  repeated  below: 


G)  is  the  harmonic  oscillatory  frequency 

L  is  the  reference  length 

V  is  the  velocity  of  imdisturbed  flow 
thenEq  7.1  becomes: 


f-1 

UJ 


+  K  -  Q(p) 


q  =  0 


(7.3) 


Eq  7.3  is  the  so-called  /^•method  equation.  It  is  the  desired  equation  for  flutter  analysis  since  its 
solution  can  provide  the  true  damping  of  the  aeroelastic  system.  However,  because  most  of  the 
imsteady  aerodynamic  methods,  including  ZONA6,  ZONA7,  ZTAIC  and  ZONA7U,  are 
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fonnulation  in  the  reduced  frequency  domain  {k  domain)  by  assuming  simple  harmonic  motion, 
the  unsteady  aerodynamic  forces  generated  by  these  methods  are  strictly  valid  only  at  zero¬ 
damping  condition.  This  implies  that  using  the  reduced-frequency  domain  unsteady 
aerodynamic  forces  (A^domain  Q{ik))  in  the  flutter  equation,  the  solution  is  valid  only  at  the 
stability  boimdary  (the  on-set  of  flutter).  Therefore,  for  reliable  damping  prediction  in  the 
complete  velocity  range  of  interest,  approximation  technique  for  flutter  solution  are  required. 

7.1  TheiT-Method 

The  basic  equation  for  flutter  analysis  by  the  AT-method  is  expressed  in  Eq  2.22,  repeated  below: 

M  +  (1  +  igJK  -  Q(/^)]q  =  0  (7.4) 


Eq  7.4  is  obtained  by  replacing  p  by  ik  in  Eq  7.3,  where  /gs  is  the  added  artificial  complex 
structural  damping  that  is  proportional  to  the  stiffness. 


The  introduction  of  /gs  was  first  proposed  by  Theodorsen  (Ref  7.1)  for  the  purpose  of  sustaining 
the  assumed  harmonic  motion.  Since  the  dynamic  pressure  can  be  written  as: 

(7.5) 


where  p  is  the  air  density. 


The  AT-method  equation  can  be  obtained  by  substituting  Eq  7.5  into  7.4  and  dividing  the  resultant 
equation  by  : 


M 


2 


Q(ik)  -  AK 


q  =  0 


(7.6) 


where: 


A  = 


0  +  ^gs) 


is  the  complex  eigenvalue  of  Eq  7.6. 

If  rigid  body  modes  exist,  Eq  7.6  cannot  be  solved  directly  since  it  contains  some  trivial  solutions 
associated  with  the  rigid  body  modes.  Therefore,  it  is  necessary  to  eliminate  these  trivial 
solutions  by  partitioning  Eq  7.6  into  rigid  body  and  elastic  modes: 


0 


Q. 

Q.r 


0  0 
0  K, 


=  0 


(7.7) 
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where  the  subscripts  r  and  /  denote  the  rigid  body  modes  and  the  elastic  modes,  respectively. 
Since: 

{q,}  =  -M„-'  M,  {q,} 

Eq  7.7  can  be  reduced  to: 

[[-M,  M„-'  M,  +  M.  ]  -  A[K„]]{q,)  =  0  (7.8) 

where: 

M.  =  M.  +  ipflj  Q, 

To  solve  for  the  eigenvalue  2.  ,  it  is  required  to  perform  the  unsteady  aerodynamic  computations 
at  several  given  reduced  frequencies.  These  reduced  frequencies  are  defined  here  as  the 
“reduced  frequency  list”.  Q(/^)  are  generated  at  a  given  Mach  number  of  interest  for  each 
reduced  frequency.  For  a  given  air  density  p ,  the  eigenvalue  of  Eq  7.6  in  terms  of  X ’s  are 
solved  in  the  complete  reduced  frequency  list.  For  n  structural  modes,  there  are  n  eigenvalues 
corresponding  to  n  modes  at  each  reduced  frequency.  The  flutter  frequency  cOf ,  the  airspeed  , 
and  artificial  damping  g*  are  given  by: 

1 

Of  = 

g,  =  Of^  Im(>^)  (7.9) 


Fig  7.1  depicts  a  typical  AT-method  results  of  the  AGARD  445.6  wing  (Ref  7.2)  with  4  structural 
modes  and  15  set  of  Q(/A:)  ranging  from  k  =  0.001  to  A:  =  0.5  computed  by  ZONA6.  The  K- 


VRe(;i) 
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method  is  performed  at  a  given  Mach  number  and  presented  in  terms  of  velocity  vs.  frequency 
diagram  (V-f  diagram)  and  velocity  vs.  damping  diagram  (V-g  diagram). 


7.1  AGARD  445.6  JT-Method  Flutter  Results  (ZONA6),  M  =  0.9,  p  =  0.000193  slug/tf 

The  V-g  diagram  shows  that  the  damping  of  mode  2  crosses  the  zero  damping  axis  at  Vt  = 
1000/sec  indicating  a  flutter  boundary  of  the  wing  structure.  The  V-f  diagram  of  Fig  7.1  clearly 
depicts  the  numerical  procedure  of  the  AT-method.  At  each  reduced  frequency,  the  flutter 
frequencies  of  all  structural  modes  are  located  along  a  radial  line  that  starts  from  the  origin.  The 
AT-method  solves  the  eigenvalues  of  the  flutter  equation,  usually  from  the  highest  to  the  lowest 
frequencies  in  the  reduced  frequency  list.  This  will  give  their  corresponding  flutter  velocities 
generally  from  the  low  speeds  to  high  speeds. 

Since  the  AT-method’s  numerical  procedure  requires  only  a  straightforward  complex  eigenvalue 
analysis  of  each  reduced  frequency,  its  solution  technique  is  efficient  and  robust.  However, 
several  drawbacks  discussed  below  make  the  /T-method  a  less  attractive  method  for  flutter 
analysis. 

•  The  solution  is  valid  only  at  g*  =0.  Other  non-zero  damping  values  are  artificial  and  may  not 
have  significant  physical  meaning. 

•  The  frequencies  and  velocities  are  computed  at  a  given  pair  of  Mach  number  and  air  density. 
This  implies  that  the  flutter  boimdaiy  computed  by  the  AT-method  generally  is  not  a  “matched 
point”  solution  in  that  the  flutter  velocity,  V^  ^  M  a„  .  The  matched  point  solution  can  be 
achieved  only  by  performing  the  flutter  analysis  at  various  air  densities  iteratively  until  the 
condition  of  Ff  =  M  a„  is  satisfied. 

•  Sometimes  the  frequency  and  damping  values  “loop”  around  themselves  and  yield  multi¬ 
value  frequency  and  damping  as  a  function  of  velocity.  This  gives  difficulty  in  tracting  the 
eigenvalue  in  the  reduced  frequency  list. 

•  The  term  —  involved  in  Eq  7.6  indicates  that  the  AT-method  cannot  generate  flutter  solution 

k 

at  it  =  0.  This  is  the  reason  why  the  AT-method  excludes  the  rigid  body  modes  from  its  flutter 
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equation.  The  failure  at  A:  =  0  also  implies  that  the  AT-method  cannot  predict  the  divergence 
speed  instability;  an  important  aeroelastic  instability  problem. 

7.2  The  P-i:  Method 

Since  its  applicability  for  flutter  method  analysis  was  first  proposed  by  Irwin  and  Guyett  (Ref 
7.3)  in  1965,  the  P-K  method  has  been  widely  adopted  by  aeroelasticians  as  the  primary  tool  for 
finding  flutter  solutions.  Hassig  (Ref  7.4)  has  given  a  detailed  description  of  the  superiority  of 
the  P-K  method  over  the  AT-method.  In  Ref  7.4,  the  equation  of  the  P-K  method  reads: 


M/  +  K  - 


{q}=0 


(7.10) 


For  simplicity,  Eq  7.10  excludes  the  structural  modal  damping  matrix,  but  it  can  be  easily 
included,  p  is  the  non-dimensional  Laplace  parameter  and  can  be  expressed  as: 

p  =  g  +  ik  (7.11) 


where; 

g  =  yk 


Y  is  the  transient  decay  rate  coefficient. 


The  P-K  method  is  an  approximate  method  of  finding  a  rate-of-decay  of  type  solution.  It  is  a 
mathematically  inconsistent  formulation  since  the  eigenvalues  p  is  expressed  as  damped 
sinusoidal  motion  while  the  Q(/^)  term  is  obtained  based  on  the  undamped  simple  harmonic 
motion.  However,  it  is  generally  believed  that  the  P-K  method  gives  a  good  approximation  of 
the  /7-method. 

Rodden  (Ref  7.5)  modified  Hassig’ s  P-K  method  equation  by  adding  an  aerodynamic  damping 
matrix  into  Eq  7. 10.  The  modified  P-K  method  equation  reads: 


Mp^ 


.K-lpy^ 


{q}  =  0 


(7.12) 


where: 

and  Q*  are  the  real  part  and  imaginaiy  part  of  Q(z^),  i.e.: 

Q(ik)  =  Q^  +  zQ'  (7.13) 

By  substituting  p  =  g  +  ik  into  the  third  term  ofEq  7. 12,  this  equation  can  be  rewritten  as: 
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(7.14) 


M/  +  K 


{q}  =  0 


By  comparing  Eq  7. 14  to  Eq  7. 10,  it  is  clearly  seen  that  the  extra  term: 


-  -  pV^^g 

2  ^  k  ^ 


is  the  added  aerodynamic  damping  matrix.  Eq  7.12  is  solved  at  several  given  values  of  V  and  p, 
for  complex  roots  p  associated  with  modes  of  interest.  This  is  accomplished  by  an  iterative 
procedure  that  matches  the  reduced  frequency  k  to  the  imaginary  part  of  p  for  every  structural 
mode.  This  iterative  procedure  is  called  the  reduced  frequency  “lining-up”  procesi’  (Ref  3.3) 
which  requires  the  repeated  interpolation  of  Q(/^)  from  these  of  the  reduced  frequency  list. 

Fig  7.2  shows  the  P-K  solutions  of  the  AGARD  445.6  wing  described  in  Sec  7.1  at  six  given 
velocities  from  700  to  1200  ft/sec  with  intervals  of  100  ft/sec.  The  damping  and  frequency 
values  at  each  given  velocity  are  presented  in  the  V-g  and  V-f  diagrams  in  Fig  7.2. 


Fig  7.2  AGARD  445.6  P-K  Method  Flutter  Results  (ZONA6),  M  =  0.9,  p=  0.000193  slug/ft^ 

Comparing  Fig  7.2  to  Fig  7.1,  it  can  be  seen  that  the  P-K  method  yields  “well-behaved”  damping 
and  frequency  curves;  no  “loop”  problem  occurs.  The  predicted  Vf  and  cOf  at  zero  damping 
agree  very  well  with  these  of  the  j^-method.  The  well-behaved  damping  and  frequency  curves 
are  believed  to  be  more  realistic  than  these  of  the  AT-method.  Also,  for  the  present  AGARD 
445.6  case,  the  P-K  method  predicts  a  flutter  mode  associated  with  the  first  mode  as  opposite  to 
the  second  mode  by  the  AT-method. 

A  more  meaningful  contrast  of  the  tow  methods  can  be  seen  when  the  divergence  speed 
instability  occurs.  Fig  7.3  shows  the  V-g  and  diagrams  of  the  /T-method  and  the  V-g  and  V-f 
diagrams  of  a  jet  transport  wing  (Ref  7.6)  with  10  modes  at  M=  0.  Both  methods  predict  a 
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flutter  speed  Fj  at  1050  ft/sec.  This  is  expected  since  at  zero  damping  both  methods  reduce  to 
the  same  equation. 

A  significant  difference  does  occur  when  a  divergence  speed  instability  appears  at  1651  ft/sec. 
This  divergence  speed  instability  is  evident  by  its  associated  zero  frequency.  The  V-g  curve  of 
the  bending  mode  predicted  by  the  iT-method  approaches  the  zero  damping  axis  perpendicularly 
but  doe  not  cross  it.  The  corresponding  damping  cure  of  the  P-K  method  has  a  discontinuity  but 
it  develops  the  divergence  speed  instability.  It  is  believed  that  this  discontinuity  is  associated 
with  the  occurrence  of  an  aerodynamic  lag  root.  In  the  next  section,  we  will  discuss  the  physical 
significance  of  the  aerodynamic  lag  root. 


Fig  7.3  Jet  Transport  Wing  at  M=0.0  at  Sea  Level  using:  (a)  K-Method  (b)  Method 
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Another  principal  advantage  of  the  cure  of  the  P-K  method  is  that  it  produces  results  directly  for 
given  values  of  velocity  and  air  density  pairs.  This  allows  the  cure  of  the  P-K  method  to  provide 
the  matched  point  solution  in  that  the  given  velocity  and  air  density  satisfy  =  M  . 
However,  in  the  next  section,  we  will  show  that  the  added  aerodynamic  damping  matrix  in  Eq 
7.12  is  vdid  only  at  small  k  or  for  linearly  varying  Q(/^).  For  highly  nonlinear  Q(/A:),  the  P-K 
method  may  produce  unrealistic  root. 

7.3  The  g-Method 

By  utilizing  a  damping  perturbation  method,  Chen  (Ref  7.7)  suggested  that  a  first  order  damping 
term  can  be  included  in  the  flutter  equation.  This  first  order  damping  term  is  rigorously  derived 
from  the  Laplace  domain  aerodynamic,  leading  to  a  flutter  solution  technique  called  the  “g- 
method”. 

Formulation  of  the  2-Method 

The  basic  assumption  of  the  g-method  lies  in  the  existence  of  an  analytic  function  of 
Q(p)  =  Q(g  +  /it),  where  g  =  yk  so  that  Q(p)  can  be  expanded  along  the  imaginary  axis 
(i.e.,  g  =  0)  for  small  g  by  means  of  a  damping  perturbation  method: 

Q(p) » Q(ii)  +  ,  forg«l  (7,15) 

% 

The  second  term  on  the  right  hand  side  of  Eq  7.15  is  generally  not  available  from  the  ^-domain 
unsteady  aerodynamic  methods.  However,  due  to  analytic  continuation,  it  can  be  replaced  by: 


gQ(p)  ^  gQ(p) 

Eq  7.16  is  valid  in  the  complete  p-domain  except  along  the  negative  real  axis  in  subsonic  flow 
(Ref  7.8).  This  implies  that  Q’{ik)  can  be  computed  from  Q(/^)  by  a  central  differencing 
scheme,  except  at  it  =  0.  At  A:  =  0,  a  forward  differencing  scheme  is  employed  to  accommodate 
the  discontinuity  of  Q(/A:)  along  the  negative  real  axis.  Substituting  Eq  7.16  into  Eq  7.15  yields 
the  approximated  /7-domain  solution  of  Q(p)  in  terms  of  k  and  for  small  g: 

Q(/7)«Q(/A:)  +  gQW  (7.17) 


fir=0 


dQ{ik) 

d{ik) 


=  Q'(/^) 


(7.16) 


Substituting  Eq  7.17  into  Eq  7.3  yields  the  g-method  equation: 


fy‘i\ 


M/j'  +  K  -  QWg  -\pV^  Qiik) 


{q}  =  0 


(7.18) 
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At  g  =  0,  both  the  g-method  and  the  P-K  method  reduce  to  the  same  form.  This  indicates  that 
both  methods  will  provide  the  same  flutter  boundary  for  zero  damping.  For  non-zero  g, 
comparing  Eq  7.18  to  Eq  7.14,  it  can  be  seen  that  the  difference  between  die  P-K  method 
equation  and  the  g-method  equation  lies  in  the  terms  Q^k  in  Eq  7.14  and  Q'iik)  in  Eq  7.18.  In 
fact,  QV^  is  a  special  case  of  This  can  be  shown  as  follows: 

Expanding  Qiik)  about  ik-Oby  Taylor’s  expansion  gives: 

Q(ii)  =  Q(0)  +  ik  Q’(0)  +  !(/*)=  Q'(0)  +  . . .  (7. 19) 

Since  all  Q'(0)  are  real,  Q(/k)  can  be  split  into  the  real  and  imaginary  parts.  It  reads: 

Qiik)  =  +  /Q'  (7.20) 


where: 

Q"  =  Q(0)  -  i*'  Q'(0)  +  ... 


(7.21) 


and: 

Q'  =  kQ'iO)  --k^  Q"iQ)  +  ...  (7.22) 

6 

Dividing  Eq  7.22  by  k  gives  the  term  qV^  in  Eq  7. 14  as: 

^  =  Q’(0)  -  4*'  Q"(0)  +  -  (’-25) 

k  6 

Differentiating  Eq  7.19  with  respect  to  ik  gives  the  term  Q'(^^)  in  Eq  7.18  as: 

Q'(ik)  =  Q'(0)  +  /itQ'(O)  +  ...  (7.26) 

Comparison  of  Eq  7,23  with  Eq  7.24  shows  that  the  equality  of  QV^  and  Q'(ik)  exists  only  if 
Q(ik)  is  a  linear  function  of  ^  or  at  ^  =  0.  This  proves  that  the  added  aerodynamic  damping 
matrix  in  Eq  7. 14  is  valid  only  if  one  of  the  above  conditions  is  satisfied.  In  fact,  if  Q(/^)  is 
highly  nonlinear,  the  P-K  method  may  produce  unrealistic  roots  due  to  the  error  from  the 
differences  between  Eq  7.23  and  Eq  7.24. 

Solution  Algorithm  of  the  ^-Method 

Substituting  p  =  g  +  ik  into  Eq  7. 18  yields  a  second-order  linear  system  in  terms  of  g: 

[g^A  +  gB  +  C]{q}  =  0  (7.25) 
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where: 


/"T/V 


A  = 


\Lj 
B  =  2/)t 


M 


C  =  -k‘ 


K 


M  -  l.pV^Q\ik)  + 


-1 

kL) 


1  fv^ 

M  +  K  -  -pV^  Qiik)  +  ik\y 

2  \L) 


Here,  Eq  7.25  is  formally  called  the  g-Method  equation.  For  completeness,  in  Eq  7.25,  we  have 
included  a  modal  structural  damping  matrix  Z.  The  solutions  of  Eq  7.25  exist  when  Im(g)  =  0. 
To  search  for  this  condition,  we  first  rewrite  Eq  7.25  in  a  state  space  form: 


[d  -  gl]{X}  =  0 


where: 


D  = 


0 

-A*C 


I 

-A‘B 


(7.26) 


(7.27) 


Next,  a  reduced-frequency-sweep  technique  is  introduced.  This  technique  searches  for  the 
condition  Im(g)  =  0  and  solves  for  the  eigenvalues  of  D  in  term  of  g;  starting  from  A:  =  0, 
incrementally  increasing  k  by  M,  and  ending  at  (^x  is  the  hipest  value  in  the  reduced 
frequency  list  of  the  imsteady  aerodynamic  computations).  The  reduced  frequency-sweep 
technique  searches  for  the  sign  change  of  the  imaginary  part  of  the  eigenvalues  between  k  and  \k 
+  Ait].  If  this  occurs,  tiie  condition  of  Im(g)  =  0  can  be  obtained  by  a  linear  interpolation  in  k  for 
the  appropriate  frequency  range.  Then  the  flutter  frequency  and  damping  2/ are  computed 
by: 


(7.27) 

k 

(7.28) 

For  0,  an  alternative  form  of  Eq  7.28  is  used  (Ref  7.9): 


Re(g) 


2r  = 


ln(2) 


(7.29) 


One  of  the  issues  in  performing  the  reduced  frequency-sweep  technique  is  the  eigenvalue 
tracking  from  kto\k+  AAj.  In  order  to  monitor  the  sign  change  of  eigenvalues,  it  is  required  that 
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the  eigenvalues  are  lined  up  at  each  k  and  \k  +  A^.  Using  the  regular  sorting  scheme  by 
comparing  the  differences  of  the  eigenvalues  at  [A  +  AAj  to  those  at  k  is  certainly  not  robust  and 
requires  small  M  values  that  may  be  costly.  This  eigenvalue  tracking  issue  can  be  resolved  by 
means  of  a  predictor-corrector  scheme. 

Predictor-Corrector  Scheme  for  Eieenvalue  Tracking 

The  predictor  predicts  the  eigenvalues  at  |A:  +  by  a  linear  extrapolation  from  the  eigenvalues 
and  their  derivatives  at  k 


g^ik  +  M)  =  gik)  +  Ak^ 

dk 


(7.30) 


da 

where  gp  is  defined  as  the  predicted  eigenvalue.  —  can  be  obtained  by  using  the  orthogonality 

dk 

property  of  the  left  and  right  eigenvectors  of  Eq  7.26.  This  leads  to; 


dk 


Y^X 


(7.31) 


where  Y  and  X  are  the  left  and  right  eigenvectors  of  Eq  7.26,  respectively,  and: 


dD 

dk 


0 

_i  dC 


0 

A-^ 


(7.32) 


Once  gp  is  given  by  the  predictor,  gp  is  used  as  the  baseline  eigenvalues  for  sorting  the  computed 
eigenvalues  at  1^  +  M|,  defined  as  gc.  The  maximum  norm  of  the  error  between  gp  and  gc  for  all 
eigenvalues  is  also  computed.  If  it  exceeds  a  certain  level,  the  predictor  could  potentially 
introduce  incorrect  eigenvalue  tracking  results  due  to  rapid  changes  of  the  eigenvalues.  In  this 
case,  the  corrector  is  activated. 

The  corrector  reduces  the  size  of  Ak  by  a  factor,  for  instance  100,  and  recomputes  gp  and  gc  at 

k  +  -^1  as  well  as  the  maximum  norm  of  the  error.  This  process  repeats  until  the  maximum 
V  lOOj 

norm  of  the  error  is  below  a  certain  level.  However,  numerical  experience  shows  that  when  the 
corrector  is  activated,  this  condition  can  be  satisfied  by  reducing  the  size  of  Ak  only  once. 
Therefore,  the  corrector  normally  would  not  increase  the  computational  time  significantly.  It 
serves  only  as  a  fail-safe  scheme. 
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At  k 


dg 


X,  —  is  also  used  to  search  for  the  condition  Im(g)  =  0  at  ^  >  k^ax  by  a  linear 
dk 

extrapolation.  Thus,  the  reduced-frequency-sweep  technique  offers  a  scheme  that  could  find  all 
real  roots  of  Eq  7.25  in  the  complete  reduced  frequency  domain. 


At  this  point,  the  issue  of  the  number  of  real  roots  that  could  exist  in  Eq  7.25  is  discussed.  For  n 
structural  modes,  the  P-K  method  and  AT-method  normally  provide  only  n  roots  of  the  flutter 
equation.  However,  as  indicated  by  Ref  26,  the  number  of  roots  could  exceed  the  number  of  the 
structural  modes  if  the  aerodynamic  lag  roots  appear.  For  instance,  if  the  exact  Theodorsen 
function  is  used,  the  number  of  aerodynamic  lag  roots  that  would  appear  is  infinite.  As  one  can 
see,  unlike  the  P-K  and  AT-methods,  the  reduced  frequency-sweep  technique  employed  by  the 
present  g-method  potentially  gives  an  unlimited  number  of  roots.  The  inclusion  of  all  activated 
aerodynamic  lag  roots  could  provide  important  physical  interpretations  of  the  flutter  solution. 
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Fig  7.4  Generalized  Aerodynamic  Forces  vs.  Reduced  Frequency  of  the  15°  Sweptback  Wing 

at  Af  =  0.45, 4  Modes 

Comparison  Between  the  s-Method  and  the  P-K  Method 


•  The  15-Degree  Sweptback  Wing  at  M  =  0. 45 

This  test  case  is  denoted  as  HA145E  in  Ref  3.3.  Four  structural  modes  are  used  for  flutter 
analysis.  The  imaginary  parts  of  the  4x4  generalized  aerodynamic  forces  matrix  (denoted  as 
Qj,)  vs.  k  are  presented  in  Fig  7.4.  Since  Im(Qy)  are  all  nearly  linear  that  gives  a  close  equality 
of  Eq  7.23  and  Eq  7.25,  the  agreement  between  the  damping  computed  by  the  P-K  method  and 
the  g-method  is  expected.  Fig  7.5  shows  the  damping  vs.  velocity  diagram  (F-g  diagram)  and 
the  flutter  frequency  vs.  velocity  diagram  (F-/  diagram)  computed  by  both  methods.  Good 
agreement  between  these  methods  is  obtained  except  the  g-method  predicts  one  extra 
aerodynamic  lag  root  (represented  by  the  crosses  in  Fig  2).  This  aerodynamic  lag  root  becomes 
active  at  F=  550  ft/sec  with  stable  damping  but  its  frequency  remains  zero.  Since  the  number  of 
roots  computed  by  the  P-K  method  is  restricted  to  be  the  same  as  the  number  of  the  structural 
modes,  at  F  =  600  ft/sec  the  P-K  method’s  reduced  frequency  “lining-up”  process  skips  the 
bending  mode  and  obtains  the  aerodynamic  lag  root;  this  creates  a  discontinuity  of  the  damping 
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associated  with  the  bending  mode  in  the  v-g  diagram  (represented  by  the  opened  triangles).  By 
contrast,  the  g-method  gives  a  continuous  damping  curve  of  the  bending  mode  (represented  by 
the  solid  triangles)  and  a  discontinuity  in  the  imping  curve  of  the  aerodynamic  lag  root  (the 
crosses)  at  1^=  550  ft/sec. 


In  order  to  investigate  how  the  g-method  obtains  the  aerodynamic  lag  root,  the  search  histoiy  in 
terms  of  eigenvalues  vs.  k  for  the  reduced  frequency-sweep  technique  is  presented  in  Fig  7.6  for 
V=  500  ft/sec  and  Fig  7.7  for  V=  650  ft/sec.  Since  there  are  4  structural  modes,  the  state  space 
form  of  Eq  7.26  provides  8  eigenvalues.  At  500  ft/sec  the  imaginary  parts  (Im(g))  of  these  8 
eigenvalues  provide  four  zero  crossings  (marked  by  the  opened  circles  in  Fig  7.6).  These  four 
zero  crossings  represent  the  four  roots  of  the  four  structural  modes. 

It  is  noted  that  the  zero  crossing  of  the  first  eigenvalue  is  obtained  by  extrapolation  from  the 
eigenvalue  and  its  derivative  at  ^  =  Amax,-  At  F  =  650  ft/sec  Im(g)  of  the  seventh  eigenvalue 
becomes  zero  at  A:  =  0  which  corresponds  to  the  occurrence  of  the  aerodynamic  lag  root.  This 
can  be  seen  clearly  in  the  expanded  view  of  Im(g)  at  small  k  (at  the  upper  right  comer  of  Fig 
7.7).  The  real  part  of  this  eigenvalue  (Re(g))  at  ^  =  0  has  a  negative  value  (Fig  7.7)  that  indicates 
this  aerodynamic  lag  root  is  stable;  however,  the  expanded  view  shows  a  potential  coupling 
between  the  aerodynamic  lag  root  and  the  sixth  eigenvalue  since  the  zero  crossing  of  the  sixth 
eigenvalue  already  occurs  at  small  k.  This  indicates  an  instability  may  appear  at  a  higher 
velocity. 
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Fig  7.5  F-g  and  F/ Diagrams  of  tihie  15®  Sweptback  Wing  at  M  =  0.45 
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Reduced  Frequency  (k) 


(a) 


Reduced  Frequency  (k) 

(b) 

Fig  7.6  Search  History  of  flie  Reduced  Frequency-Sweep  Technique  at  500  ft/sec,  (a) 

Imaginary  Damping  and  (b)  Real  Damping 
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Fig  7.7  Search  History  of  the  Reduced  Frequency-Sweep  Technique  at  V-  600  ft/sec 


•  Jet  Transport  Wing  at  M=  0.0  with  10  Modes 

This  test  case  has  been  discussed  in  Sec  7.2.  Ten  stractural  modes  are  used  for  flutter  analysis 
but  only  the  results  of  the  first  bending  and  torsion  modes  are  presented  in  the  V-g  and  V-f 
diagrams  shown  in  Fig  7.8.  Two  types  of  instability  are  predicted  by  both  the  P-K  method  and 
the  g-method:  a  flutter  speed  at  V=  1056  ft/.sec  and  a  divergence  speed  at  V-  1651  ft/.sec.  This 
agreement  is  expected  since  at  g  =  0  the  flutter  equation  of  both  methods  reduce  to  the  same 
form.  Three  aerodynamic  lag  roots  are  obtained  by  the  g-method  and  their  frequencies  are  all 
zero  throughout  the  velocity  range  of  interest.  Both  of  the  first  and  second  aerodynamic  lag 
roots  become  active  at  the  same  speed  (V  =  1400  fl/sec).  After  this  speed,  the  second 
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aerodynamic  lag  root  forms  a  super-stable  mode  but  at  this  speed  the  damping  of  the  first 
aerodynamic  lag  root  jvunps  suddenly  firom  zero  to  -  0.1  then  gradually  crosses  the  zero-damping 
axis,  forming  a  divergence  type  of  instability  at  V=  1651  fl/sec.  At  this  divergence  speed,  the 
third  aerodynamic  lag  root  becomes  active  and  suddenly  jumps  to  a  high  value  of  unstable 
damping  (Fig  7.8).  This  is  an  interesting  phenomenon  because  it  indicates  that  this  divergence 
speed  could  be  a  bifurcation  point.  Determining  the  third  aerodynamic  lag  root  is  bifurcated 
from  the  first  aerodynamic  lag  root  or  originates  on  its  own  needs  further  investigation. 

Similarly  to  the  first  test  case,  the  damping  curve  of  the  bending  mode  computed  by  the  F-K 
method  has  a  discontinuity  while  that  of  the  g-method  remains  a  smooth  curve.  The  damping 
curve  of  the  torsion  mode  computed  by  both  methods  are  in  excellent  agreement.  The  frequency 
curves  of  the  two  structural  modes  computed  by  both  methods  also  in  good  agreement  except  for 
the  absence  of  the  three  aerodynamic  lag  roots  of  the  F-K  method. 
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Fig  7.8  V-g  and  ^-/Diagrams  of  the  BAH  Wing,  Af  =  0.0, 10  Modes 
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•  Two  Degrees  of  Freedom  Airfoil  at  M=  0.0 

This  test  case  is  adopted  from  Ref  7.10  and  is  derived  from  the  case  denoted  as  HAMS  A  in  Ref 
3.3  but  with  the  fuselage  grid  point  being  constrained.  The  center  of  gravity  is  located  at  37% 
chord.  Fig  7.9  presents  the  variations  of  the  2  x  2  Q,;  vs.  k.  In  this  case.  Fig  7.9  shows  that  the 
imaginary  parts  of  Q,y  is  not  linear.  Therefore,  some  difference  in  flutter  results  between  the  P-K 
method  and  the  g-method  is  expected.  First,  for  clarity,  the  V-g  diagram  computed  by  the  g- 
method  alone  is  presented  in  Fig  7.10.  Two  aerodynamic  lag  roots  are  found.  Again,  it  seems 
that  the  second  aerodynamic  lag  root  is  bifurcated  from  the  first  one  at  F=  210  ft/sec  where  a 
divergence  instability  occurs.  The  comparison  of  the  damping  and  flutter  frequencies  between 
the  P-K  method  and  the  g-method  is  shown  in  Fig  7.11;  however,  for  clarity,  the  second 
aerodynamic  lag  root  is  not  repeatedly  shown.  In  Fig  7.1 1  the  results  computed  by  the  transient 
method  are  also  presented.  The  transient  method  is  based  on  a  time-domain  unsteady 
aerodynamic  method,  therefore  it  can  be  considered  as  a  p-method.  All  of  the  three  methods 
predict  the  same  instabilities:  a  divergence  instability  at  F=  210  ft/sec  and  a  flutter  instability  at 
V  =  250  ft/sec.  The  damping  ciuves  of  the  first  and  second  modes  computed  by  the  g-metiiod 
correlate  well  with  those  of  the  transient  method.  But,  again,  the  P-K  method  gives  a 
discontinuous  damping  curve  of  the  first  mode. 

For  the  case  of  the  center  of  gravity  moved  to  45%  chord,  the  F-g  diagram  shown  in  Fig  12.a 
indicates  that  the  flutter  instability  (at  F=  170  ft/sec)  occurs  before  the  divergence  instability  (at 
V—225  ft/sec).  Again,  this  is  well  predicted  by  all  three  methods.  The  frequency  curves  in  the  V-f 
diagram  (Fig  12.b)  computed  by  the  g-method  show  a  similar  trend  as  those  of  the  transient 
method.  But  the  curves  of  the  P-K  method  are  discontinuous  at  F  =  100  ft/sec  where  an 
aerodynamic  lag  root  appears  (not  obtained  by  the  transient  method  but  well  captured  by  the  g- 
method).  This  results  a  poor  correlation  of  the  V-f  curves  obtained  by  the  P-K  method  with  the 
other  two  methods. 


Fig  7.9  Generalized  Forces  of  2  D.O.F.  AirfoO,  C.G.  @  37%  Chord  (HA145Al),3/=  0.0, 2  Modes 


134 


Fig  7.10  2  D.O.F.  Airfoil,  C.G.  @  37%  Chord  (HA145Al),ii/=  0.0, 2  Modes 
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Fig  7.11  Damping  and  Frequency  vs.  Velocity  of  2  D.O.F.  Airfoil,  C.G.  @  37%  Chord  (HA145A1), 

M  =  0.0, 2  Modes 
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Fig  7.12  Damping  and  Frequency  vs.  Velocity  of  2  D.O.F.  Airfoil,  C.G.  @  45%  Chord  (HA145A2), 

Af  =  0.0, 2  Modes 


•  Three  Degrees  of  Freedom  Airfoil  at  M=  0.0 

This  test  case  is  denoted  as  HA145A  in  Ref  3.3.  A  fuselage  free-free  plunge  mode  is  added  in 
the  above  two  degrees  of  freedom  case.  The  V-g  and  V-f  diagrams  for  the  case  of  the  center  of 
gravity  located  at  37%  chord  are  shown  in  Fig  7. 13  and  those  for  45%  chord  are  in  Fig  7. 14.  For 
both  cases,  the  so-called  “dynamic  divergence”  (Ref  7.11)  occurs  and  its  speeds  and  frequencies 
are  well  predicted  by  all  three  methods:  the  P-K  method,  the  g-method,  and  the  transient  method. 
Both  the  g-method  and  the  transient  method  capture  one  aerod3mamic  lag  root  (in  the  45%  chord 
case,  the  g-method  obtains  a  second  lag  root  but  it  becomes  active  at  the  dynamic  divergence 
speed  and  is  not  discussed  here).  Unlike  the  restrained  structures  of  all  previous  test  cases  where 
the  frequency  of  the  lag  roots  remains  zero,  the  aerodynamic  lag  root  of  the  present  unrestrained 
structure  ‘takes  off  from  the  zero-frequency  axis  then  couples  with  the  bending  mode.  This 
coupling  of  the  lag  root  and  bending  mode  forms  a  “dynamic  divergence”  instability. 
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As  indicated  by  Ref  7.11,  this  dynamic  divergence  has  a  non-zero  frequency  and  could  be 
defined  as  a  low-frequency  flutter  instability.  On  the  other  hand,  the  P~K  method  generated  lag 
root  somehow  refuses  to  ‘take  off  from  the  zero-frequency  axis.  This  problem  of  the  P-K 
method  is  probably  due  to  the  fact  that  since  Qj,  of  the  present  test  case  is  nonlinear,  the  P~K 
method  is  valid  only  at  ^  =  0  for  non-zero  damping.  This  k-Q  condition  restricts  the  frequency 
of  the  lag  root  from  being  a  non-zero  value  and  results  in  a  poor  correlation  in  the  V-f  diagram 
with  the  other  two  methods. 
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Fig  7.13  Damping  and  Frequency  vs.  Velocity  of  3  D.O.F.  Airfoil,  C.G.  @  37%  Chord  (HA14SA2), 

A/’=0.0,3Modes 
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Fig  7.14  Damping  and  Frequency  vs.  Velocity  of  3  D.O.F.  Airfoil,  C.G.  @  45%  Chord  (HA145A2), 
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Fig  7.15  Johnson  Configuration  Generalized  Aerodynamic  Forces,  M=0.84, 17  Modes 


138 


•  The  Johnson  Configuration  at  M=  0.84  with  1 7  Modes 

This  test  case  is  adopted  from  Ref  7.12.  The  Johnson  configuration  has  three  rigid  body  modes 
and  14  elastic  modes.  The  imaginaiy  parts  of  Q,y  vs.  k  for  i  and  y  =  4,  5  and  6  presented  in  Fig 
7.15  show  that  spikes  occur  at  small  k.  The  cause  of  the  spikes  is  probably  due  to  poor 
aerodynamic  modeling;  but  this  is  not  an  issue  to  be  discussed  here.  Since  Q,y  are  highly 
nonlinear,  a  large  difference  between  the  results  obtained  from  the  P-K  method  and  the  g-method 
is  anticipated.  In  fact,  in  this  case  the  P-K  method  breaks  down  (Ref  7.12)  and  its  results  are 
totally  unreliable.  It  is  believed  that  the  break-down  of  the  P-K  method  is  caused  by  the 
unrealistic  roots  produced  by  the  nonlinear  Q,y.  In  order  to  validate  the  g-method  result,  the  K- 
method  is  used  for  comparison. 
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Fig  7.16  Damping  and  Frequency  vs.  Velocity  of  Johnson  Configuration,  M  =  0.84, 17  Modes 

There  are  13  aerodynamic  lag  roots  obtained  by  the  g-method.  Due  to  the  spikes  at  small  k, 
some  of  them  become  active  even  at  very  low  speed.  These  lag  roots  are  not  presented  here.  Fig 
7.16  shows  the  V-g  and  V-f  diagrams  obtained  by  the  AT-method  and  the  g-method  for  the  first 
three  elastic  modes;  denoted  as  mode  4,  5,  and  6.  It  can  be  seen  that  both  methods  predict  the 


same  flutter  boundary  around  V  =  470  ft/sec.  The  good  agreement  between  the  AT-method  and 
the  g-method  indicates  the  robustness  of  the  g-method’s  solution  algorithm. 


It  is  generally  believed  that  the  AT-method  is  only  valid  at  the  g  =  0  condition.  The  present  work 
also  proves  that  the  P-K  method  is  valid  at  the  conditions  of  g  =  0,  ^  =  0,  or  =  0  ,  where 


n  >  1.  The  g-method  generalizes  the  AT-method  and  the  P-K  method.  It  is  valid  for  all  k  and  up  to 
the  first  order  of  g.  This  first  order  term  of  g  is  rigorously  derived  from  Q(p)  by  a  damping 
perturbation  method. 


The  present  work  also  provides  a  theoretical  foundation  for  the  g-method  that  can  be  used  to 
estimate  the  error  of  large  damping  (beyond  the  first  order  assumption)  due  to  the  truncation  of 
the  higher  order  terms  of  g.  However,  based  on  the  formulation  of  the  g-method,  adding  higher 
order  terms  in  g  seems  to  be  straightforward. 
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